NAVAL POSTGRADUATE SCHOOL 
Monterey , California 




THESIS 

l 1 



DYNAMIC STALL ANALYSIS UTILIZING 
INTERACTIVE COMPUTER GRAPHICS 

by 

Eric L. Pagenkopf 
* ° * 

March 19S8 



Thesis Advisor 
Co-Advisor 



MX. Plat/er 
LAV. Carr 



Approved for public release; distribution is unlimited. 



T 23 91 22 



I’nclassified 

sec j r it> classification of this page 



1 3 a T\pe of Report 


13b Time Cc\ered 


IJ Date of Report ivmr. month, day) 


1 5 Page Count 


Master's Thesis 


From To 


March 1988 


90 



REPORT OOf.LMLM AIION R \(,I 



la Rrp^n Secuuv. Classification T ncla>>ihed 



2a Scour U> Chwdieation Atubopp 



2b aao n P o \\ \vz i r/i * n g S ched u le 

i P-miming Ogpruziumn Reror Nurrbens) 



*■> 1 \a "j <4 roi torn-iii— CM:. n mon 

Na\ al Postgraduate School 



Cl Ik 1 ; Symbol 

( if afvlkable > 67 



\d j.ess f r./rv. Z/P . ae ) 

Montorev. C.\ 0^945.51100 



Sa Name of Funding Sponsoring Organiz 



anon 



6b Office Symbol 
t if Gpph' able > 



6c Address y city, state. ana ZIP coaej 



\b Rev.iv \c Mrr.Vnus 



3 Distribution \ wtwiibdU; A Report 

Approved lor public release: distribution i? uiduniwd. 



5 Monitor Orgar\;v r, '’n R.-pm’. Nunccmi 



“a Nan : uf M. nu’rmg Ore iiuzanoi: 
Naval PoMgiaduute School 



7h Add: ess , :,Ty. ym.v cvia 1 ZIP cc^e [ 
Montercv. CA '4VM3oijnii 



9 Procurement instiumeiu Identification Nurnt 



nber 



in ^ ,i jrce oi f undine >’« 



>mcer i 



Pr^gr.in Element \o | Project N r | r ? k \ 0 | \\ . ■ l •* w 



11 Title ( include security d^i/c^tion; DVNAM1C S TALL ANALYSIS UTILIZING IN TERACTIVE COMPL I FRGRAPIL 
ICS 



12 feisonal Authors) Liic L. Pagenkopf 



in Supplementary Notation The views expressed m tliis thesis are those of the author and do not reilect the official policy or po- 
sition of the Department of Defense or the L’.S. Government. 



1 " Cosatt Cedes 


Field 


Group 


Subgroup 















IS Subject Terms 1 continue on re\ene ij necessary and identi/v by block number) 
dynamic stall. computer graphics IRIS.llow visualization 



lv Abstract i continue on reverse if necessary and identify by blxk number) 

A Navier-Stokes pioblem solver, developed by L. N. Saiikar, is modified to provide dynamic, interactive giaphical pres- 
entations of predicted flow field solutions about a NACA-00I2 airfoil section, oscillating in pitch, m order 10 demonstrate the 
capabilities of dynamic graphics applications in the study of complex, unsteady flows. Flow held solutions in the loirn oi 
pressure coefficient and stream function contour plots about an airfoil experiencing dynamic stall me plotted utilizing an [PCS 
3000-series woikstation and Graphical Animation System (GAS) software, developed by Sterling Software for NASA. These 
full cycle solutions, in conjunction with dynamic surface pressure disiiibution plots and integrated hit. pitching moment and 
drag coefficient data, are compared to existing experimental data in Older to provide an indication of the validity ot the code s 
iur-field solution. Full procedural documentation is maintained in order to pro\ide an efficient analysis tool for me 111 future 
oscillating airfoil studies planned by the NASA-Ames Fluid Mechanics Laboratory and the Naval Postgraduate Schooi De- 
partment of Aeronautics and Astronautics. 



20 Distribution Availability of Vos tract 
53 unclassified unlimited D sarr.3 as report 



□ one users 



22a Name of Responsible Individual 
M.F. Platzer 



21 Abstract Security Classification 

Unclassified 



22b Teleohone include Area code) 

(40S) 646-23 1 1 



22c Office Symbol 

54Ss 



DD FORM 1473.84 MAR 



S3 APR edition may be used until exhausted 
All other editions aie obsolete 



securitv classification of f h:>' race 



Uild issilied 



1 



Approved for public release; distribution is unlimited. 



Dynamic Stall Analyse tilizing 
Interactive Computer uraphics 

by 

Eric L. Pagenkopf 
Lieutenant. United States Na\y 
13. A., California State University, Northridge, 1980 

Submitted in partial fulfillment of the 
requirements for the degree of 

MASTER OF SCIENCE IN AERONAUTICAL ENGINEERING 

from the 



NAVAL POSTGRADUATE SCHOOL 
March 1988 



ABSTRACT 



A Navier-Stokes problem solver, developed by L. N. Sankar, is modified to provide 
dynamic, interactive graphical presentations of predicted flow field solutions about a 
NACA-0012 airfoil section, oscillating in pitch, in order to demonstrate the capabilities 
of dynamic graphics applications in the study of complex, unsteady flows. Flow field 
solutions in the form of pressure coefficient and stream function contour plots about an 
airfoil experiencing dynamic stall are plotted utilizing an IRIS 3000-series workstation 
and Graphical Animation System (GAS) software, developed by Sterling Software for 
NASA. These full cycle solutions, in conjunction with dynamic surface pressure dis- 
tribution plots and integrated lift, pitching moment and drag coefficient data, are com- 
pared to existing experimental data in order to provide an indication of the validity of 
the code's far-field solution. Full procedural documentation is maintained in order to 
provide an efficient analysis tool for use in future oscillating airfoil studies planned by 
the NASA-Ames Fluid Mechanics Laboratory and the Naval Postgraduate School De- 
partment of Aeronautics and Astronautics. 
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I. INTRODUCTION 



It has been recognized for some time (Ref. 1] that an airfoil oscillating in pitch to 
angles of attack greater than the static stall angle will surpass the traditional stail bunier 
and generate normal forces which exceed those attainable in the static case, lias dy- 
namic stall phenomenon is attributed to the aft movement of a strong shed voitex along 
the upper surface of the aiifoil, carrying with it an induced velocity held which radically 
and dynamically changes ehordwise pressure disti ibutions. In geneial, an accurate in- 
terpretation of the dynamic stall mechanism will significantly impact a variety of appli- 
cations, all of which involve dynamic lifting suifacc motions and unsteady flew 
separation. Originally, dv naniic stall analysis efforts were directed toward helicopter 
aerodynamics, where sharp increases in oscillatory torsional loading and thus, blade 
stress, can reduce the fatigue life of rotor mechanical components and vortex-induced 
aerodynamic loading can generate adversely phased pitching moments, resulting in stall 
flutter [Ref. 2). Much current interest concerns the feasibility of exploiting dynamic stall 
forces for effective, sustained maneuvering in the high angle of attack. "post-staU" flight 
regime, or supermaneuverability, for next-geneiation fighter and attack aircraft [Ref. 3). 
Historically, attempts to analyze such complex, unsteady behavior relied heavily or. 
empirical data, obtained from often laborious and time-consuming tests. With the ad- 
vent of the supercomputer, however, computation of actual viscous flow fields about 
moderately complex computational models can now be numerically achieved in a matter 
of minutes through solution of the Reynolds-averaged Navier-Stokcs equations [R.ef. 
4|. These solutions, in and of themselves, however, are not sufficient to promote insight 
into the mechanics and physics involved in such flows. This additional requiiement. for 
effective visual portrayal of the flow field, is satisfied by application of high performance 
interactive computer graphics workstations and associated software. 

The long range goal of the Computational Fluid Dynamics field is the development 
of a thoroughly verified computer code for unsteady aerodynamics which, among a 
myriad of other applications, will provide future aircraft designers with the opportunity 
to derive full advantage from application of the dynamic stall phenomenon. To this end. 
the Fluid Mechanics Laboratory (F.ML) of the XASA-Amcs Research Center (ARC), 
in conjunction with the Naval Postgraduate School Department of Aeronautics and 
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Astronautics has planned an assortment of parallel and complementary oscillating airfoil 
windtunnel experiments and computer simulations. 

The current study utilizes existing dynamic graphics packages to generate real-time 
projections of flow fields about a NACA-0012 airfoil oscillating in pitch and experienc- 
ing deep dynamic stall, as provided by a Navier-Stokes solver developed by L. N. Sankar 
[Refs. 5.6]. A modified version of the Sankar code was submitted via the FML front end 
VAX, to the NASA-ARC Cray X-.V1P 48 computer, which output flow field solutions 
at specified intervals throughout the oscillatory cycle. From this data, graphics files were 
generated which, in turn, were submitted to the Graphics Animation System (GAS) 
software as developed by Sterling Software under contract to NASA-ARC, and run on 
an IRIS 3000-series graphics workstation. (Final output, for demonstrative purposes, 
was then transferred to video tape.) Thus, the results of the study were two-fold. First, 
procedures by which interactive computer graphics could be efficiently (in terms of both 
computer-time and man-hours) and effectively (in terms of information display) incor- 
porated as an analysis and verification tool for future FML studies, were developed, 
tested and refined. Secondly, the resultant graphics were utilized as a tool for the on- 
going verification of the Sankar code. 
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II. DESCRIPTION OF THE SANKAR CODE 



Simulation of complex phenomena occurring in real fluid Hows requires accurate 
solutions to the full Navier-Stokes equations. The Sankar Navier -Stokes solver, devel- 
oped for Blade-Vortex Interaction (BVh studies, solves the two-dimensional, unsteady, 
compressible Euler and Navier-Stokes equations in strong conservation form, utilizing 
an alternating direction, implicit method as the time marching algorithm. A body-fitted 
C-grid. with clustering in the normal diiection is utilized to discretize the flow field. 
Turbulent shear stresses are simulated with a two-layer algebraic eddy-viscosity model, 
with modifications as described later. The code may thus be utilized for solution of 
steady or unsteady, inviscid or viscous and laminar or turbulent flows. 

A. GOVERNING EQUATIONS 

In Cartesian coordinates, the 2-D, unsteady, compressible Navier-Stokes equations 
in strong conservative form may be written as 



d t q + S X E + OyF = Re '(t 5 X R -T SyS) (2.1) 

where q, E, F , R and S are four element vectors with entries corresponding, in order, to 
the equations of continuity, momentum (x- and v-) and energy. If non-dimensionalized 
such that all values of length, velocity (u and v), density ( p ) and total energy per unit 
volume (e) are normalized with respect to section chord length (c), free stream speed of 
sound ( a x ), free stream density ( p x ) and p x a x1 respectively, these vectors may be writ- 
ten as 
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The Reynolds number (Re), pressure (p), speed of sound (a) and stress terms ( O are 
defined as follows. 
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Re = 



(2.3) 




P = (y - !)[<?- 0 .5p(u 2 + v 2 )] 



(2.4) 



a 2 = y(y - 1 )0/p - 0.5(w 2 + y 2 )] 



(2.5) 



t xx = (/ + 2u)u x 4- / Vy 



(2.6) 



T 



= /4 w >. + v a) 



r yy ~ (^' 4* 2 J u) v y + 

I 2 

P 4 = z/r ;cc + ut jcv 4- &/V (y — l)5 x a 

5 4 = wt x> , + vr^ + kPr~\y — l)<5yz 2 



Assuming Stokes' hypothesis to be valid, the constant ). is defined as -2 3 n . Since only 
airflow is considered, the ratio of specific heats ( y ) is defined as 1.4 and the Prandtl 
number (Pr), as one. 

Neglecting viscous terms results in the Euler equations, as the right-hand side of 
equation (2.1) goes to zero. L’sc of the strong conservative form (i.e., the continuity 
entry in £, S(pu)ldx , is solved in its present form, vice the more mathematically correct 
form of pSu/Sx + uSp/dx) allows the identical conservation of physical quantities (mass, 
momentum and energy) when finite difference schemes are applied. 

B. TRANSFORMED GOVERNING EQUATIONS 

The flow field region must be discretized into a transformed, finite difference mesh, 
or computational plane, in order to allow numerical solution of the governing equations. 
The most efficient grids are rectangular in shape with regular spacing or spacing gradi- 
ents. They must, however, correspond to a flow field grid which provides high resolution 
(minimal spacing) in regions where gradients are large, particularly in the boundary layer 
and about the leading edge. In response to the coordinate transformation involved in 
the grid generation process, the governing equations, too, must be transformed. By de- 
fining ( and rj as functions of the cartesian coordinates (time is not transformed), this is 
simply a 2-D mapping procedure accomplished by application of the transformation 
Jacobian, 
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the governing equation becomes 



O.q + < 5 ? £ + < 5 / = Re~\d ; R + 5, t S) 



where 



q = q\J 



E=(Z t q +Z x E + t; y F)IJ 



E= (n,q + n x E + n F)lJ 



R = {(, X R +Z y S)IJ 

S = (q x R +VyS)jJ 



T XX - V- + 2/t)(C x W ? + V x U n ) + /(C,V’{ + Vy v n ) 

*xy = + >1yU n ) + (C^ + q x v n )~] 

T yy = (;. + 2^)(C y r c + ^) + /(C x m + r] x u n ) 

Ri = + vr xy + kPr~\y - lXC^a 2 + q x S n a 2 ) 

= UT xy + viyy + kPr~\y - l)(C/5 c a 2 + >l y S n a 2 ) 



(2.7) 



( 2 . 8 ) 



(2.9) 



( 2 . 10 ) 



( 2 . 11 ) 
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C. GRID GENERATION 

Grid generation for discretization of the How Held region may be accomplished 
conformal mapping, algebraic methods or numeiically solving a set of partial differential 
equations. The original Sankar code utilizes either of the latter two methods, while the 
version currently vectorized for use in this study utilizes only the the algebraic method, 
which results in a body-fitted, sheared parabolic coordinate system or C-grid. Utilization 
of such body-fitted grids allows synchionous rotation of the entire grid and aii foil sec- 
tion for dynamic cases, using simple tiigonometric relations. This coordinate sestem 
satisfies the general gi id requirements for smoothness throughout and line spacing in 
regions where high gradients exist, such as the boundary layer or leading edge. 

Normalized geometric airfoil shape data, in Cartesian coordinates, is input in table 
format to the code, which utilizes an interpolative procedure to compute additional 
points and smooth the surface. The trailing edge region is modelled as a \ortex sheet 
shape, or "cut", which smoothly leaves the airfoil, tangent to the mean camber line at 
the trailing edge. Algebraic manipulation of the section surface and cut allows grid 
generation in the transformed c - ;/ plane. Figure 1 shows the resultant grid when 
mapped back to Cartesian coordinates. 

Uniform spacing in the transformed plane's c -direction results in fine spacing at the 
leading edge in the real plane, but coarse spacing in the trailing edge region, due to 
uniform increments across the axis (as opposed to the standard C-grid which results in 
a region of finer spacing at the trailing edge). ;/ -spacing in both planes increases, ap- 
proximately exponentially, in directions normal to the airfoil surface, with the magnitude 
of initial spacing being user defined. Though easy to generate, this grid's effectiveness 
in the analysis of certain flow cases has proven somewhat limited, due to its coarseness 
in the trailing edge region [Ref. 6, p.44|. 

D. APPLICATION OF BOUNDARY CONDITIONS 

Consideration of an airfoil impulsively started from rest in a fluid with uniform 
properties throughout, provides the initial conditions for solution of the parabolic Euler 
and Navier-Stokes equations (Eq. 2.9) in the computational plane. In addition, bound- 
ary conditions and artificial dissipation terms are required in order to achieve accurate 
solutions. 

Three boundaries exist which are of interest, the surface, the Tar-field boundary (grid 
limit) and the trailing edge cut. Boundary conditions at the surface are driven by the 
no-slip condition which, in response to viscous effects, requires the fluid velocity at the 
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Figure 1. Algebraic C-Grid in Cartesian Coordinates 
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boundary to equal the boundary velocity. Thus, in this reference frame, u and v are set 
to zero on the solid surface. Since the surface is considered adiabatic, both b<:,'tb/ and 
oc/S'l are set to zero and. likewise, pressure distributions are determined from 6pi0>; and 
oplo'i equal to zero at the solid boundary. (This is equivalent to neglecting stress con- 
tributions in the momentum equations, which is acceptable when dealing with high 
Reynolds numbers flows.) Boundary conditions at the cut are driven by the necessity 
to avoid discontinuities in the "continuous" poition of the How field. Since the grid is 
extremely dense in the >i -direction in this region, averages of the values of the two 
nearest interior points are assigned to points along the cut, allowing a smooth transition. 
Since the grid cannot economically be generated large enough to avoid disturbance \e- 
locities at the fur-field boundary, boundary conditions must account for the presence of 
the airfoil. Linear small disturbance theory is applied to determine perturbation veloci- 
ties at points along the boundary, which are then added to free stream conditions. 
Downstream boundaries are treated such that entropy changes can be convected out of 
the computational domain [Ref. 6. p. 29), allowing shocks and boundary-layer generated 
vorticity to pass through the grid. 

It has been found that spatial derivatives are sensitive to the decoupling of odd and 
even points which necessarily occurs in central difference schemes. This results in the 
generation of high frequency errors in regions of large pressure gradients, such as are 
present in shocks or about stagnation points. Due to the high Reynolds numbers en- 
countered in these flows, dissipation provided by the viscous terms are not enough to 
eliminate the errors, necessitating the addition of two artificial dissipation terms, em- 
bedded in the numerical schemes. An explicit artificial viscosity term is input by the 
user, with a proportional implicit term assigned by the code. 

E. TURBULENCE MODELLING 

Since during the derivation of the Navier-Stokes equations, no assumptions are 
made regarding flow type, these equations are instantaneously valid for both laminar 
and turbulent flows. The large range of time and spatial scales encountered in turbulent 
flows, however, makes solution of the instantaneous Navier-Stokes equations viitually 
impossible, at this time. As a result, the equations are time- or ensemble-averaged. Use 
of such equations, in response to turbulence, results in additional turbulent shear stress 
terms or the Reynolds stresses (named for Oswalde Reynolds, who initiated turbulence 
studies in the 1880 s), which effectively are increases in shear stresses due to turbulent 
motion. The difficulties associated with the existence of additional unknowns without 
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any additional equations is described as the closure problem and is dealt with through 
the use of turbulence models. These are based on empirical knowledge of turbulent 
Hows and thus, provided solutions will be approximate and require some confirmation 
from experiment [Ref. 7]. 

The most common tuibulence modelling method involves mixing lengths, as devel- 
oped by Prandtl. Assuming that turbulent fluctuations are essentially the result of ve- 
locity perturbations between adjacent streamlines, and that all fluctuating velocity 
components at a given point are of the same magnitude, it can be shown that turbulent 
or eddy viscosity ( f.t T ) is proportional to the magnitude of the local vorticity ( cu ) [Ref 
S, p. 3S7|. In Prandtl's mixing length model, the proportional term is the square of a 
characteristic length, related to fluid turbulence intensity. Prandtl suggested this char- 
acteristic length ( / ) be treated as 



/ = ky 



where y is the normal distance from the fluid boundary and k is empirically obtained. 
Thus, the key to obtaining accurate solutions when utilizing models of this type lies in 
the mixing length expression. 

The Baldwin- Lomax two-layer eddy viscosity model, based on Cebeci's two-layer 
model, is presently used in the Sankar code. This model divides the boundary layer into 
two regions, an inner layer and an outer layer, with separate methods for determining 
eddy viscosity used in each. The boundary’ for the two layers is defined as that point 
where eddy viscosities produced by the two methods match. 1 he inner layer utilizes a 
mixing length model where eddy viscosity starts from zero at the wall and is defined by 
the following. 

(«rU = P' ! l<»l (2.12) 

The mixing length term is defined by 



1 = k>-D (2.13) 

where y is the normal distance from the wall, k is the Von Karman constant and D is the 
Van Driest damping function, given by 

D = [l - exp( VMP] (2.14) 



where A~ is an empirical constant. The outer layer eddy viscosity model is defined by 
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outer FC C pp ^wake^Kleb^y^ (-.15) 

where K and C ep are constants. The K motf intermittencv function, given by 

F Kleb O’) = [l + 5'~(CKteby'ym ax) 1 (--1^) 



ensures that eddy viscosity approaches zero as the edge of the boundary layer is ap- 
proached and the flow assumes external characteristics. C KUb is a constant. F„ akt is a 
function defined the following relation. 



F\vake nu n (V m a x f max’ ^wk} 



max t-diflF max) 



(2.17) 



Uj, f is the magnitude of the velocity profile's velocity range. F mAX is the maximum value 
provided by the following. 



F(y)=y\uj\D (2.18) 

y max is the value of y corresponding to F miX . Constants in the Sankar code are defined 
as follows: 

A + = 26.0 



k = 0.4 



AT = 0.016S 

Qp -1.6 
C-Kleb = 0.3 

c wk = 0.25 

In order to account for turbulence outside the boundary layer, such as that which 
occurs in the dynamic stall process, a modified turbulence model is available which ef- 
fectively increases the outer model's mixing length. In this case, F max and y max are deter- 
mined by the following. 

F{y)=y 2 \o>\D (2.19) 



10 



The modified turbulence model was developed in response to early predictions of flow 
separation at high angles of attack and has been found to cause premature reattachment 
during the downstroke of dynamic cycles. Its use, therefore, is only recommended until 
stall onset [Ref. S.p. 33]. 

F. CODE STRUCTURE 

The code is structured such that the following user defined parameters are input 
from logical unit 05 (appended to the end of the code for Cray processing): grid size, step 
size (dt >, artificial viscosity magnitude, mean angle of attack (a 0 ), oscillation magnitude 
(aj. angle for suspension of the modified turbulence model, reduced frequency (k). free 
stream Mach number (M x ), Reynolds number (Re), distance of the first >/ - contour from 
the airfoil wall, starting time, pitch and restart flags and airfoil geometry data. Added 
to this list for the present study are the number of time steps for the code to march on 
the present run, the number of steps between each plot (output interval) and the total 
number of steps and plots completed on all previous runs. Variables are then initialized, 
previous stored solution datasets are retrieved and read (to allow incorporation of all 
datasets, including those from the present run, into a single, combined file) and flow field 
starting conditions are input. A series of subroutine calls then generate the grid 
(AIRFOL, SING, TABINT, WRAP), cluster grid points for viscous flows (CLUSTR. 
STRTCH), rotate the airfoil (ROTGRID) to its actual angle ol'attack and compute the 
initial metrics (METRIC). 

At this point, the code is fully initialized for commencement of the flow field sol- 
ution process, which is conducted via an iterative loop. At each iteration, time is first 
marched forward one time step, followed by computation of the time dependent values 
of angular velocity, angle of attack and step change of angle of attack, according to the 
relations: 



to = 2kXf 00 sin(2A.V/ 00 r) 



( 2 . 20 ) 



a,- = «(,-«, cosilkM^t) 



oq_i = «o - a, cos[ 2 kMJyt - dt)~\ 
da. = a — oq_j 

The grid is then rotated, in the physical plane, by applying the the following relations 
at each grid point. 
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x = x cos(da) — y sin(</a) 



( 2 . 21 ) 



y = y co s(tfo.) — x sin (dsavepha) 

Following recomputation of the metrics, the solution is computed by subroutine SLPS, 
the ADI-algorithm, which calls DISSIP, for computation of the explicit dissipation 
terms, STRESS and RESI, for computation of the inviscid and viscous terms, respec- 
tively, AM ATI and MATRIX l, for computation of the Jacobian and its inverse in the 
C -direction and AMAT2 and MATRIX2, for computation of the Jacobian and its in- 
verse in the rj -direction. Subroutine STRESS also calls subroutine EDDY, the turbu- 
lence model, for computation of the viscosity coefficient. WALLBC enforces the wall 
boundary conditions and finally, solution files, as discribed later, are generated. 

Prior to resumption of the loop, performance coefficients are generated by subrou- 
tines LOAD and CPPLOT. Surface pressure coefficients are obtained from the relation: 

= P^-Poo — ( 2 . 22 ) 

1 /2p oe ( A/ooO 

where p b is the surface pressure and p x is the free stream pressure. Skin friction coeffi- 
cients are a function of the wall shear stress (rj according to the relations: 







r w =Re~ i MJ V 



(2.23) 



\V=U y - V x S J(M n Xr + y ? ) 

The aerodynamic loads of lift, drag and moment about the quarter-chord, are computed 
by the following relations. 



Q = C n cos(a) — C f sin(a) (2. 24) 

C d = C„ sin(a) — C r cos(a) 

C m = I C p L{y- y C j A )dy + (x - x c/4 )dx] 
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C, and C, are the normal and tangential forces obtained from summing surface pressure 
and skin friction forces according to the following. 



G. CODE MODIFICATIONS 

The following modifications were made to the code, resident on the FML VAX. 
during the course of this study. 

1. Solution output commands are placed within the loop in order to provide interval 
output. 

2. Read commands are placed at the beginning of the main program in ol der to allow 
storage of all solutions (previous restarts and the current run) in a single combined 



3. The airfoil section is rotated to the initial angle of attack, vice rotation of the free- 
stream direction. 

4. Restarts may be made from Plot3D (output) format as well as vector format. 

5. Additional Plot3D files (surface pressure and skin friction line plots) are output 
from subroutine CPPLOT. 

6. Grid dimensions are made true variables. 

7. User inputs are reformatted for ease of use. 





file. 
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III. DYNAMIC GRAPHICS GENERATION 



A. INTERACTIVE COMPUTER GRAPHICS 

Due to the complexities and time dependent phenomenon associated with unsteady 
flows, attainment of acceptable clarity in flow field solutions requires complete descrip- 
tive information across fine grids, at a large number of minute time intervals. Thus, 
analysis and visualization of data generated by codes such as Sankar's, is difficult due 
to both the volume of information provided and its temporal nature. Requirements 
identified for effective flow field visualization include maximization of the bandwidth of 
information transfer, such that it closely matches the capabilities of' the human eye, 
maximization of the quality of graphical displays, by enhancing key features and sup- 
pressing others, and maximization of the controllability of information [Ref. 9j. Addi- 
tional advantages in information transfer can be achieved through the use of redundant 
coding and structured displays [Ref. 10). In response to these requirements, interactive 
computer graphics workstations have been evolved to complement the super-computer. 
Workstation capabilities, in terms of geometrical transformation and screen update dis- 
patch have been utilized by programmers to produce effective representations of flow 
field motion, through svnehronization of coordinated data sets. Solution claritv is en- 
hanced by the high degree of spatial resolution and large range of colors afforded by the 
workstation displays. The interactive capabilities which currently exist for semi-complex 
flows serve to improve visual cues for display of three-dimensional data sets and allow 
immediate access to regions of interest. 

B. HARDWARE 

The NAS A- ARC Fluid Mechanics Laboratory is currently equipped with an 
IRIS-3000 series workstation configured with 4Mbytes of display list memory and 
equipped with z-clipping and z-buffering hardware and a multiple key mouse. The IRIS 
display processing unit's refresh memory is arranged as a two-dimensional array (768 x 
1024), with row-column positions matching display pixel x-y coordinates. 24 bits are 
reserved for each pixel at eight bits (256 intensities) per color, resulting in a range of 
16,777,216 possible colors. When displaying dynamic graphics, the refresh memory' is 
divided, in order to meet user demands, since flyback time (about 1.3 n sec) is too short 
for complete picture update. This double buffering mode utilizes half the memory for 
refreshing the screen displays, while the other half is dynamically updated. This ensures 
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smooth motion on the display, but divides the number of bitplanes per pixel, in half. 
Thus, the number of colors available drops to 4096. In order to avoid a similar decrease 
in the range of colors available, color maps are utilized. In this case, pixel values in the 
refresh memory are routed to an index or color look-up table (with each entry composed 
of 24 predefined bits, which define the color), vice direct routing to the intensity digital- 
to-analog converter. The IRIS transformation rate, via a single, composed transforma- 
tion matrix in homogeneous coordinates, is S0,o00 coordinates per second, with points 
and lines generated at 3 million pixels per second, or 40,000 inches per second. Polygons 
can be filled (flat shading) at the rate of 44 million pixels per second, or 70,000 square 
inches per second. These capabilities allow generation of most displays at the rate of 
approximately 10 screens per second, which, under ordinary conditions, is greater than 
the fusion frequency of the human eye. Interactive transformations are ordinarily ac- 
complished via the mouse. 

C. PLOT3D SOFTWARE 

Solutions provided by the Sankar code consist of multiple binary files, specifically 
formatted for input into the Plot3D graphics software, authored by Pieter Buning for 
NASA. Two filetypes exist, which are labeled as XYZ or grid tiles and Q or How 
quantity files. Both types are three-dimensional matrices, with placement of data within 
the matrices corresponding to the row-column addresses of mesh points. XYZ-files have 
a depth of two in the two-dimensional case and provide cartesian coordinate definitions 
of each mesh point. Q-files have a depth of four in the two-dimensional case and pro- 
vide computed flow quantities at each mesh point, corresponding to the q-vector of 
equation 2.1. Headers for each file list free-stream Mach, Reynolds number, angle of 
attack and time. 

Additional flow' field quantities at each mesh point are computed according to the 
following relations. [Ref. 11] 

Definitions: 



y — 1.4 



R = 1 .0 (Gas constant) 



M = Vic 
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Density: 



P = ?U) 

Poo * 1-0 



(3.1) 



Pressure: 



Temperature: 



Enthalpy: 



Energy: 



P o 






1 + 



-1 



Ar 




P = (>' - OpC^o - -5E 2 ] 



1 





T _ P p 
loo 



r ° = 7 0 + i_ r~‘' /2 ] 



h = y[e 0 - -5V 2 1 



K = y^ooo 

, ^ p 

h 0 = e 0 + ~ 




(3.2) 



(3.3) 



(3.4) 



(3.5) 
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Entropy: 



= 



<7 * T) 
P 



' oo 

e 0oo ^/so ^ ^ oo 

^ = ,5E 2 



5 = 111 



[ 



P 

Poo 




= 0 



Pressure Coefficient: 






P ~ P ec 

■SPooK 




Po — Poc 



■ 5 P oo V 'oo 



1 



(3.6) 



(3.7) 



Based on these values, the Plot3D software will produce plots of scalar functions 
(density, pressure, temperature, enthalpy, energy, velocity, entropy, momentum and 
shocks) suitable for contour plots, vector functions (velocity, vorticity, momentum and 
pressure gradient), particle trace functions (particle traces and vortex lines), shock wave 
locations and grid functions (computational meshes and walls). 

Particle traces are generated using trilinear interpolation of values of the vector 
function inside a computational cell and second-order Runge-Kutta steps to advance the 
particle in space. Five steps are required for each cell. The algorithm for computing 
shocks computes the Mach number component in the direction of the local pressure 
gradient. Locations where this value decreases through 1.0 is plotted as a shock. Two- 
dimensional stream functions are calculated by integrating the mass flow across a coor- 
dinate line. 



17 



Output from Plot3D is available in a variety of formats and is a function of user- 
defined attributes. For dynamic plots, output is in the form of graphics files, formatted 
for further manipulation by animation software. Device independent plotting (DIF) is 
enabled through use of the ARCGraph binary file format which allows data to be ma- 
nipulated by a collection of function libraries and utilities developed by the Advanced 
Computer Research Center at NASA-ARC. Parameter data, along with graphics prim- 
itives (device independent op-codes) are contained in these files, allowing data manipu- 
lation by the Cray, IRIS and attached VAXes. Plot3D generation of graphics files 
requires that several attributes be defined. In order to reduce user-tasking, Plot3D will 
initially search for a filename-specific initialization communications file, which contains 
all required inputs. Input files must be read singularly, which makes necessary' the sep- 
aration of those combined plotting files received from the Cray. Output graphics files, 
however, are once again stored in combined files in order to speed further processing. 
User-defined attributes include axis scales and extreme values, function, number of 
contours, line colors and types and wall definitions. 

Subroutine CPPLOT, within the Sankar code, provides grid (XYZ) files which con- 
tain plotting information for surface pressure coefficient and skin friction coefficient line 
plots, corresponding to each output flow field solution file. These line plots are scaled 
and translated after separation from the combined files and plotted utilizing the grid 
function. Likewise, an angle of attack pointer plot, utilizing a sine wave format (gener- 
ated external to the code), is also provided. Q-files containing dummy variables (not 
utilized for grid plot generation) are also provided to complement the line plot 
XYZ-files, as required by Plot3D. 

D. GRAPHICS ANIMATION SYSTEM SOFTWARE 

The final step in the visualization process is the animation or synchronization of 
coordinated graphics datasets. The Graphics Animation System (GAS) is a software 
package developed by Sterling Software under contract to NASA and provided to edu- 
cational institutions throughout the country. It is device specific, running only on IRIS 
workstations, and is written in the C programming language under the UNIX operating 
system. Graphics files, in ARCGraph format, are read and stored in the IRIS' display 
list memory as objects by the menu-driven program and assembled according to se- 
quence file instructions (stored transformation matrices). In order to smooth motion 
associated with geometric transformations, up to 30 linearly interpolated matrices are 
provided by the program, between user-identified keyframes. Titles, legends and object 
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attributes are also available. Ultimately, flow Held representations may be tumsleived 
to video tape or 16mm film. Interactive viewing (ordinarily available in real-time) is 
controlled via a mouse, within the menu's "view data window. 

A complete description of GAS and associated video hardware may be found in 
Reference 12. 
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IV. RESULTS AND DISCUSSION 



A. PROCEDURAL CONSIDERATIONS 

Integration of the FML IRIS with related CFD software and the Cray X-MP/48 
was tested by plotting the flow field about an airfoil experiencing deep dynamic stall, as 
provided by the Sankar code. The advantages of plotting dynamic stall are two-fold. 
First, it is an extremely complex and time-dependent phenomenon which is difficult to 
visualize from discrete data sets. Thus, the accuracy of solutions for this t\pe flow field 
can only be determined by consideration of the complete flow field in both time and 
space. Dynamic graphics are therefore ideally suited for its study. Secondly, the FML 
is currently heavily involved in studies of dynamic stall, which include verification of the 
Sankar code. Cray output is thus used to full advantage. 

In order to provide smooth animations, an extremely large number of plotting sets, 
provided at equal time intervals, must be generated by the code during the oscillatory 
cycle. Since access and transfer of this data is required several times during the graphics 
generation process, it was found that storage of datasets in combined files provided the 
greatest efficiency, especially in transfers between the Cray and IRIS. Software resident 
on the IRIS is then employed to separate the combined file into individual datasets. 
(Embedded in these programs are schemes to scale data, as necessary to control their 
ultimate on-screen positions, allowing, for example, placement of line plots in a specific 
corner of the display or plotting, for comparison, two flow fields side-by-side.) Individ- 
ual datasets are automatically assigned names which correspond to those contained in 
Plot3D initiation files, also resident on the IRIS. These files currently contain com- 
mands for the processing of 50 datasets. Thus, the combined file separation and 
ARCgraph file production (Plot3X) process must be completed at intervals. The use of 
repetitive file names during this process, however, increases automation to such an ex- 
tent that only a limited number of user inputs are required to complete the process. 
ARCgraph files as provided by Plot3D are again in combined file format, allowing easy 
(a single user command) input into animation software. 

B. THE DYNAMIC STALL PROCESS 

A prerequisite to review of any code-provided flow field solution is consideration of 
available empirical information. The following dynamic stall characteristics have been 
observed for both harmonically oscillating airfoils [Ref. 13] and airfoils undergoing 
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monotonic angle of attack increases [Refs. 14, 15|. First, the airfoil stalls at an angle of 
attack which exceeds the static stall angle. Second, the corresponding normal forces and 
pitching moments exceed those attainable in the static case and, finally, a rapid change 
in pitching moment magnitude, termed as "moment stall" by Harris and Pruyn [Ref. 16] 
occurs several degrees in azimuth prior to the normal force decrease, or "lift stall", vice 
simultaneous occurence as in the static case. The process which results in these char- 
acteristics can be broken down into a series of chronological events, as depicted in Fig- 
ure 2 from Reference 13. While the specific characteristics for a given airfoil are a 
function of free-stream Mach number. Reynolds number, reduced frequency, mean angle 
of attack and oscillation amplitude, empirical data obtained from a NACA-0012 airfoil 
at k = .15, Re = 2.5xl0 6 and a = 15’ — 10’ sin(cor) is provided in the following discussion 
to illustrate the dynamic stall process. 

At point (a) of Figure 2, the static stall angle is exceeded without any detectable 
chanse in flow over the airfoil. The boundarv laver remains thin with no evidence of 
flow reversal at the surface. At point (b), ( a = 19 — 20’) flow reversal appears at the 
surface. Over the majority of the airfoil, the boundary layer remains thin and attached. 
At the rear portion of the airfoil, however, the boundary layer thickens as the flow 
gradually decreases to zero velocity. At point (c), large eddies appear in the boundary' 
layer. By point (d), flow reversal has spread over much of the chord, from the trailing 
edge to x.'c = 0.3. With as much as 50% of the airfoil experiencing flow reversal, no 
changes in C„ or C m are detected. At point (e), ( a = 23.4’ ) a vortex forms. The 
boundary layer at the front of the airfoil abruptly breaks down (simultaneously from x c 
= 0.0 to 0.3) and the vortex moves downstream at approximately 35 - 45% of the free- 
stream velocity. At point (0 the lift-curve slope increases beyond the 2 na limit of 
quasi-static flows. By point (g), the pressure distribution is altered sufficiently to 
produce a noticeable divergence in the pitching moment (moment stall) but is accom- 
panied by a continued increase in lift. Maximum lift occurs at approximately mid-chord, 
followed by a sharp decrease (point (h)), which identifies lift stall. (Boundary layer sep- 
aration has occurred to such an extent that continued increases in angle of attack will 
not result in continued increases in lift.) At point (i) ( a = 24.95 0 ), maximum negative 
moment occurs. The vortex moves off the trailing edge at a = 24.8 ° (downstroke) and 
C„ and C m move toward static levels. At point (j), the airfoil is fully stalled. At point (k), 
boundary layer How begins to reattach, progressing from the leading edge at approxi- 
mately 25 - 35% free-stream velocity and is complete by a = 7 0 (downstroke). At point 
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PITCHING MOMENT. C M NORMAL FORCE. C 



(h) 




(k) 




(a) STATIC STALL ANGLE EXCEEDED 



(b) FIRST APPEARANCE OF FLOW 
REVERSAL ON SURFACE 




(c) LARGE EDDIES APPEAR IN 
BOUNDARY LAYER 




(d) FLOW REVERSAL SPREADS OVER 
MUCH OF AIRFOIL CHORD 




(e) VORTEX FORMS NEAR 
LEADING EDGE 




(f) LIFT SLOPE INCREASES 




(g) MOMENT STALL OCCURS 




(h) LIFT STALL BEGINS 



(i) MAXIMUM NEGATIVE MOMENT 

(j) FULL STALL 




(k) BOUNDARY LAYER REATTACHES 
FRONT TO REAR 




(I) RETURN TO UNSTALLED VALUES 



I iguie 2. I lie Dynamic Stall Process (from Carr, Ref. 13) 
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(1), after some hysteresis ( a = 6 0 , upstroke) the separated region has fully closed and 
unstalled values of C n and C m are reestablished. 

C. PROCEDURAL VERIFICATION 

The case chosen for procedural verification closely corresponds to the Carr data of 
the previous section and exactly matches conditions previously run when the code was 
initially vectorized for use on the ARC Cray [Ref. 17]. There are several advantages 
associated with this duplication. First, it allowed storage of a complete record of this 
case (data saved at 300 intervals, vice 12) for future access. Second, it provided a simple 
means by which to determine if code modifications were correctly incorporated. Finally, 
it provided a baseline, in conjunction with empirical data, against which future sensitiv- 
ity studies could be compared. 

Inputs for this case were as follows: 

Re = 3.45x 10 6 
.\C = .2S3 
St = .005 sec. 

>] -min = .00005 
co = .151 

Artificial Viscosity = 5 
Grid size: 157 x 40 

Oscillatory Cycle: a = 15° — 10° cos(cor) 

The original Baldwin-Lomax turbulence model was used throughout the oscillatory 
cycle. Figures 3 through 7 are coefficient of pressure contour plots of the predicted How 
field. (Dynamic plots are similar to these Plot3D-provided plots, but do not contain 
explicit quantitative information.) As previously noted [Ref. 17, p. 48], under these 
conditions, moment coefficients are consistently underpredicted as compared to exper- 
imental data. Drag and lift coefficients closely match experimental data below approxi- 
mately 15 and 18 degrees angle of attack, respectively. Early prediction of flow 
separation, however, forces inaccuracies in quantitative solutions beyond these values 
and underprediction of the maximum lift coefficient obtainable. These findings highlight 
the rationale behind the modified turbulence model. Information provided by the dy- 
namic plots, shows general qualitative agreement with experimental data, during the 
upstroke. A vortex forms near the leading edge of the airfoil section and progresses 
down its upper surface, gradually increasing in size. As the downstroke begins, however, 
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Figure 3. Pressure Contours, a = 15° - 10° (.151t), a = 5.00°, 5.24 °, 5.92°, 7.02°, 
8.47°, 10.23° , (top to bottom, left to right), = .283, Re = 3.45 x 10*. 

^7 min = -00005. 
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ure 4. Pressure Contours, a = 15° - 10° (.151 1), a = 12.20°, 14.30°, 16.43°, 
18.50°, 20.40°, 22.07° , (top to bottom, left to right), M. x = . 2S3, 
Re = 3.45 x 10\ »/ mm = .00005. 
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Figure 5. Pressure Contours, a = 15° 
24.98°, 24.60°, 23.80° . (top 
Re = 3.45 x 10\ ?/ mm = .00005. 





10° (.15 It), a = 23.41°, 24.36°, 24.89°, 
to bottom, left to right), M x = .283, 
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Figure 6. Pressure Contours, 2 = 15° - 10° (.151t), 2 = 22.59°, 21.03°, 19.20°, 

17.19°, 15.07°, 12.9-4° , (top to bottom, left to right), t \f oe = .2S3, 
Re = 3.45 x 10\ ;/ mm = .00005. 
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Figure 7. 



Pressure Contours, a = 15° - 10° (.151 1), a = 10.92°, 9.07°, 7.50°, 6.27°, 
5.43°, 5.00° , (top to bottom, left to right), = .233, Re = 3.45 x 10\ 



;/ min = .00005. 
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rather than shedding from the trailing edge, the vortex stagnates at approximately the 
65 % chord point and gradually dissipates. Line plots of surface pressure coefficients 
more accurately portray (as compared to traditional carpet plots) the turbulent effects 
associated with the existence of the vortex. During the vortex dissipation phase, these 
line plots show a slight (gradually decreasing) pressure rise in the vicinity of the vortex. 
The subtle differences between these plots and those which would have been obtained 
had the vortex actually shed, are lost when this information is integrated for lift, drag 
and, to a lesser extent, moment coefficient data. Thus, quantitative data may not be 
eood indicators of flow field solution accuracy. 

D. SENSITIVITY ANALYSIS 

The following code parameters which should effect dissipation of the vortex have 
been identified. 

1 . Artificial viscosity magnitude. 

2 . Turbulence model parameters. 

3 . Grid resolution. 

A limited sensitivity analysis into the effects of grid resolution changes in the trailing 
edge region was conducted in order to provide direction for follow-on studies. During 
the initial phase of this analysis, only qualitative results will be considered, for those 
reasons mentioned above. Dynamic graphics, therefore, are ideally suited for this task, 
especially when data scaling capabilities are utilized in order to plot comparative data- 
sets, side by side. 

The grid is relatively coarse in the region in which the vortex becomes stationary and 
dissipates, as compared to those regions in which it behaves as expected. In the first 
sensitivity case, therefore, grid dimensions were held constant and the distance of the 
first constant- 17 line from the airfoil surface was increased to .001 , resulting in a finer 
grid in the trailing edge region. Figure S shows the effects of this change on the vortex. 
Within this grid, the vortex once again becomes stationary, but effectively splits, such 
that a secondary vortex is shed while the original vortex dissipates. While this does not 
correspond to any process encountered in actual flows, it does indicate that final sol- 
utions are sensitive to this parameter. 

Changes in 17 -spacing have two effects on the flow field solution. Grid spacing in 
the trailing edge region becomes more fine at the expense of grid spacing in the boundary 
layer, which becomes more coarse. The extent of influence of the boundary layer sol- 
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Figure 8. 



Pressure Contours, a = 15° - 10° (.15 It), a = 
21.85 , 19.20°, 16.13° , (top to bottom, left 
Re = 3.45 x 10\ >/ min = .001. 



24.89°, 24.85°, 23.80°, 
to right), M x = .283, 
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ution on the vortex, or far field solution, is not yet clear. It is known. ho\ve\er. that 
boundary layer solutions are extremely sensitive to >/ -spacing parameter changes, re- 
quiring an initial value of .00005 for accurate solutions (Ref. 18], Further study, there- 
fore, is required in order to determine whether vortex behavior changes observed in the 
previous case were in response to the altered boundary la\er solution (which would re- 
commend investigation of turbulence model parameters) or the altered grid. In order to 
facilitate this study, the existing grid would require enlargement (from 161 \ 41). Un- 
fortunately, this results in such a large memory allocation from the Cray, that it is pro- 
hibitively inefficient. Replacement of the entire grid generation process is therefore 
currently under consideration for this code. 

E. CONCLUSIONS 

The current study has resulted in completion of the following items. 

1. Full procedural documentation (Appendix A) has been developed for efficient, 
code-independent, two-dimensional, dynamic graphics production on FML and 
NFS IRIS workstations. Data must be in Plot3D format and may be provided by 
either code or experiment. A method for simultaneous viewing of multiple datasets 
has been incorporated. 

2. The Sankar Navier-Stokes solver (Appendix B) has been modified 10 allow dynamic 
graphical presentation of its flow field solutions. 

3. Dynamic graphics applications in the study of complex flows were illustrated by 
means of an introductory sensitivity analysis, using the Sankar code. This limited 
analysis identified either grid resolution in the trailing edge region or the boundary 
layer solution as parameters to which vortex behavior is strongly dependent. (This 
suggests that incorporation of a more sophisticated grid generation scheme would 
be worthwhile.) 

In order to derive full benefit from application of this technology, existing and future 
empirical data should be stored in plotting format, allowing the availability of visual re- 
cords of multiple flow’ functions. Also, modification of additional Navier-Stokes solvers 
for dynamic output will allow effective eode-to-code and code-to-experiment compar- 
isons, all in an effort to eventually develop an accurate flow field predictor. 



31 



APPENDIX A. PROCEDURES FOR GENERATING DYNAMIC 

GRAPHICS AT FML 



A. THE ANIMATION PROCESS 

Creation of dynamic graphics for flow Field visualization at the NASA-Ames Re- 
search Center Fluid Dynamics Laboratory involves Five basic processes: 

1. Navier-Stokes solver is submitted to the Cray X-MP, 48 and dynamic output saved 
in the Central Storage Facility (CSF) in combined Files. 

2. Combined Files are converted to IRIS binary and transFerred to an FMLIRIS1 ac- 
count. 

3. Combined Files are separated into individual plotting Files on the IRIS and con- 
verted to combined graphics (.gra) Files by Plot3D. 

4. Graphics Files are Fed to the Graphics Animation System software (GAS) and dy- 
namic graphics are plotted. 

5. Resultant Files and plots are transFerred to storage or video tape. 

Prior to beginning the animation process, the code must be modiFied to provide dy- 
namic output in Plot3D format. This involves: 

1. Placement of solution output commands within the iterative loop along with some 
interval flag. 

2. IF restarts are to be used, placement of read statements at the beginning of the 
program to allow storage of all provided solutions in a single combined File. 

3. Storage of any additional required plotting information (line plots) in Plot3D XYZ 
Files for eventual display utilizing Plot3D's grid Function. 

Full documentation of Plot3D Formats and Functions may be Found in Sterling Software 
technical note 15, A Guide to Generating Movies using Plot3D and GAS, by Greg Howe. 
The procedures outlined below refer specifically to a version of the Sankar Navier-Stokes 
solver, modified For dynamic output. They can easily, however, be generalized to apply 
to any code (or experiment) From which dynamic output is desired. While many meth- 
ods of completing these steps may exist, those outlined below have proven to be the 
most time efficient. 

1. Vax Submitted Cray Jobs 

A Full dynamic cycle, utilizing a 161 x 41 grid, requires over 5000 seconds of 
Cray time, making job chaining necessary to obtain a reasonable priority. With 900 
second jobs, 24 to 4S hours are usually sufficient For the Full computation to be com- 
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pleted. Since significant differences exist between the first JCL from which Plot3D files 
are saved and subsesquent restart JCL's, two JCL's are maintained on identical copies 
of Sankar's code. 

INITIAL. JCL. Accesses CSF or Cray tape (via STAGEX commands) for initial 
solution, which may be in either PlotiD or matrix format. Combined Plot3D files are 
initialized on CSF and assigned PDN's and ID's. Tape 05 input data (appended to end 
of code) is assigned in accordance with definitions included in MAIN portion of the 
code. (If restart is from dynamic data, AOA will not be a whole number, TSTARf 
must be adjusted accordingly, in order to have accurate time 'AOA matches.) 

RESTART. JCL. Tape 05 inputs must identically match those from 
INITIAL. JCL, with the exception of TSTART.CSTP.CPLT and possibly FORMAT. 
INITIAL. JCL-defmed PDN ID solutions are accessed (read, stored and appended with 
additional solutions). saved (with identical PDN ID s) and scrubbed (oldest versions 
deleted). The initial restart solution is accessed with library routine GETP3D. written 
by Greg Howe, which will read a combined file, skip a specified number of datasets 
(DSSK1P) and access the next dataset (or datasets, if NUMDS is not equal to one, the 
default). Thus, if 

DSSKIP + CPLT - 1, 

the correct solution will be provided for restarts. Subsequent restarts illustrate the ad- 
vantage of combined file name (PDN ID) repetition, as they simply require the following 
RESTART .JCL and Tape 05 updates: 

1. DSSKIP 

2. Job chain commands. 

3. CSTP 

4. CPLT 

Once the optimum number of plotting sets for a cycle is determined, plotting 
files for an AOA pointer on IRIS displays can be generated using Plot3D. (500 sets 
maximum, or modify dimension statements.) 

Line Plots. Any pointer or line plot information may be generated, either within 
the code or externally, and stored in XYZ files for further processing. Utilization of 
Plot3D function 0 and proper descriptors of walls will produce the desired plots. 
SINE.TXT is a program which generates plots for a dynamic AOA pointer in sine wave 
format. In this program. Tape 05 includes values for pointer scaling (axis length with 
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no scaling is In ) and positioning on the IRIS screen and initial pointer location. For 
example, if the cycle starts from the minimum AOA, 

PNLOCI = .75(NPLOTS). 

Hardcopy printouts may also be obtained lioin these grid files by use of pro- 
grams such as CPPLOT.TXT. This program provides printouts of upper and louer 
surface pressure coefficient and skin friction values at x c locations, including Cp plot 
generation. 

2 . Path to the IRIS 

SEND.JCL. Utilizes GETP3D and SEXDWKS to convert to IRIS binary and 
transfer complete or partial combined files, or individual plotting sets, from the CPvAY 
to the IRIS. (Since I O's between the Cray and IRIS can get "clogged" if too many 
transfers are requested, the individual method is not recommended.) This may need to 
be completed in stages to avoid exceeding the IRIS' memory capacity. 

3. Creating Graphics Files 

Combined files on the IRIS are separated into individual plotting sets by the 
programs GETSIN.F (also generates the dummy Q file, "qsin.iris"). GETCP.F (also 
scales and positions Cp or skin friction on the IRIS screen) and GE’I’X.F (for XYZ and 
Q flow field datasets). These programs require some initial user inputs prior to running. 
Output dataset names will automatically match Plot3E) initialization file names. Since 
the Plot3x run for each dataset type (X, P and SIX) requires specific arguments, initial- 
ization files for each data type exist separately in files XIXT.COM, PIXT.COM and 
SIXT.COM. Since Plot3X searches for the file name "PLOT3DIXT.COM ', these files 
must be renamed appropriately, utilizing the "mv" conmiand. Entering the command 
"plot3x" will then generate multiple graphics files which are stored in a combined file 
matching the initial Q file name encountered by Plot3D, (i.e., "qUUl.gra"). This file is 
then stored on an appropriate subdirectory or fed to GAS. 

4. Notes on GAS 

For movies, the maximum number of objects per sequence is fifty. The 
"seq(n).seq" files are 50 object far-field solution files (seq(n).seq = objects n(l-50)), 
which use object 4u0 for titles. They can be fed to GAS (after object input) via the AUX 
I O window. For interactive viewing, the VIEW DATA window is the only usable op- 
tion. (Xo sequence files are required in view data.) Full GAS documentation is avail- 
able in the GAS User's Guide, available from Steiliug Software. 
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5. Storage 

Since IRIS memory space is severely limited at FML. resultant graphics files 
should be stored elsewhere, if not in immediate use. Downloading to 1 4" cassette tape 
is easily accomplished (also a recommended backup), as is transfer to CSF. 
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B. SUPPORT CODES 
1. INITIAL. JCL 



JOB. JOFLYNAVY. T-900. ERIC PAGENKOPF , X4269 . 
ACCOUNT , AO . US- . UPW- . 

• , 

CFT.Ot^A.OFF-S. 



• RESTART COMmIANDS: 

CSF . ACCESS . CSFACCT- . CSFPSWO . DN-OLDS LN . PDN-NSE509 . I OVA IDES . 
ASSIGN . DN-OLDSLN , A-FT07 . 

• . ■■■'■■ 11 1 ~ 1 

•.SAVE DATA FROM PRESENT RUN IN COMBINED FILE: 

ASSIGN , DN— XYZ , A-FT20. 

ASSIGN, DN=Q, A-FT21. 

ASS I GN , DN— CPX , A-FT50. 

ASSIGN , DN— CFX , A-FT60. 

• . 

LDR . MAP-PART , SET-ZERO . 



CSF . SAVE . CSFACCT- , CSFPSWO , Df^XYZ . 
CSF . SAVE . CSFACCT- , CSFPSWO . D*=0 , 
CSF , SAVE . CSFACCT- . CSFPSWO . DN-CPX . 
CSF . SAVE , CS FACCT- . CSFPSWO . DN-CFX . 



PDN-SNKR05X, 
PDN—SNKR05Q . 
PDN-SNKR05P . 
PDN-SNKR05F , 



I OFLYNAVY . 
I OFLYNAVY . 
IOFLYNAVY. 
IOFLYNAVY . 



• .ACCESS. DN=SENDVAX,PDN=SENDVAX. I OSTTRDM.OWN-RFTRDM. 

• .SENDVAX.DN=XYZ.VDN= , RAL"JIANPS password": : [ J I ANPS . PAGAN . DATAjXT 1 .DAT' . 

• .SENDVAX.DN=Q.VDN= , RAL"JlANPS password": : [ J I ANPS . PAGAN . DATA]QT 1 .DAT’ . 

• .ACCESS, DN-SENDWKS.PDN-SENDWKS. IOSTTRDM.OWN— RFTRDM. 

• . SENDWKS . DN-03X49 . VDN— ’ FMLIRIS1 " j i aa password” : : ”/u/ j i aa/pagan/q . i r i s" ’ , NOREC. 
• .SENDWKS. DN-CFX. VDN- ’FML I RISr’j i aa password": :"/u/j i aa/pagan/cf 01 .dot"’ .NOREC. 

•.JOB CHAINING COMMANDS: 

FETCH . DN-J0B2 . TEXT- ’ [ PAGAN ]NSMULT . TXT ; 2 * . 

REWIND, DN-J0B2. 

SUBMIT. DN-J0B2. 



/EOF 
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REST ART. J CL 



JOB. JN-FLYNAVY.T-900. ERIC PAGENKOPF , X4269 . 
ACCOUNT . AO . US- . UPW- . 



CFT.ON— A.OFF-S. 



• .RESTART COMMANDS: 

• . CSP . ACCESS . CSFACCT- . CSFPSWD- . DN-OLDSLN . PDN-NSE509 , 
CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-GETP3D . PDN-GETP3D . 



CSF . ACCESS , CSFACCT- , CSFPSWD- . DN-XYZ 1 . 

CSF, ACCESS. CSFACCT-. CSFPSWD-. DN-01 . 

CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-CPX1 . 

CSF . ACCESS . CSFACCT- , CSFPSWD- . DN-CFX 1 . 
REWIND. DN-Q1 . 

GETP3D . 1—01 .OOLDSLN. DATATYPE-O. DSSK I P-49. 
REWIND. DN-OLDSLN. 

REWIND. DN-XYZ1 . 

REWIND. DN-01 . 

REWIND. DN-CPX1 . 

REWIND. DN-CFX 1 . 

ASSIGN.DN— XYZ1 , A-FT90 . 

ASSIGN, DN-01 . A— FT 9 1 . 

ASSIGN. DN-CPX1 . A-FT92 . 

ASS I GN. DN-CFX 1 . A-FT93. 

ASS IGN. DN-OLDSLN. A-FT07. 



PDN— SNKR05X . 
PDN-SNKR05Q . 
PDN-SNKR05P . 
PDN-SNKR05F . 



ID-VALDES . 
ID-RFRICC. 
ID— FLYNAVY. 
I D— F LYNAVY . 
ID— FLYNAVY. 
I D-F LYNAVY. 



• . — 1 — ■ ■■ ■ 

•.SAVE DATA FROM PRESENT RUN IN COMBINED FILE: 
ASSIGN. DN-XYZ. A— FT 20 . 

ASSIGN. DN=0. A— FT 2 1 . 

ASS I GN . DN—CPX , A— FT 50 . 

ASSIGN. DN-CFX. A— FT 60 . 



LDR , MAP-PART . SET-ZERO . 



CSF. SAVE. CSFACC T- 
CSF , SAVE , CSFACCT— 
CSF. SAVE. CSFACCT- 
CSF. SAVE. CSFACCT- 
CSF, SCRUB. CSFACCT 
CSF, SCRUB. CSFACCT 
CSF. SCRUB. CSFACCT 
CSF. SCRUB. CSFACCT 



.CSFPSWD-. DN-XYZ. PDN-SNKR05X. 
.CSFPSWD-.DN-Q. PDN-SNKR05Q, 

. CSFPSWD-, DN-CPX, PDN-SNKR05P, 
.CSFPSWD-. DN-CFX, PDN— SNKR05F, 

- . CSFPSV*)— . PDN-SNKR05F, I D-F LYNAVY . 
. CSFPSWD- . PDN— SNKR05P . I D-FLYNAVY . 
. CSFPSWD- . PDN-SNKR05X , I D-FLYNAVY. 
. CSFPSWD- . PDN-SNKR05Q . I D-FLYNAVY . 



I D-FLYNAVY. 
I D-FLYNAVY. 
I D-FLYNAVY. 
ID— FLYNAVY. 



. ACCESS . DN-SENDVAX . PDN-SENDVAX .ID-. OWN- . 

. SENDVAX , DN— XYZ . VDN- ’ RAL" J I ANPS password": : [ J I ANPS . PAGAN . DATA]XT1 .DAT’ . 

. SENDVAX . DN— Q . VDN— ' RAL" J I ANPS password": : [J 1ANPS. PAGAN . DATA] QT 1 .DAT' . 

. ACCESS . DN-SENDWKS . PDN-SENDWKS .ID-. OWN- . 

. SENDWKS , DN=^33X49 , VDN-’ FMLIRIS1 " j i oa password" : : "/u/ j i aa/pagon/q .iris”’, NOR EC . 
. SENDWKS. DN-CFX. VDN- ’FMLlRISV'j iaa password": :"/u/j i aa/pogon/cf 01 .dot'” .NOREC. 



•.JOB CHAINING COMMANDS: 

FETCH , DN— J0B3 , TEXT— ’ [PAGAN] J0B3 . TXT ' . 
REWIND. DN-J0B3 . 

SUBMIT . DN-J0B3 . 



/EOF 
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3. SEND.JCL 



JOB. JN-GETX.T-30. ERIC PAGENKOPF, X4269 . 

ACCOUNT . AO . US= , UPW- . 

• . 

CSF . ACCESS . CSFACCT* . CSFPSWO , DN=XYZ . PDOSNKR38X . I OF LYNAVY . 

CSF , ACCESS , CSFACCT* . CSFPSWO . DN*Q . PDN=SNKR38Q . I OF LYNAVY . 

• . 

CSF . ACC ESS . CSFACCT= , CSFPSWO . DN=GETP3D . PDN-GETP3D . I ORFR ICC . 

ACCESS . DN=SENDWKS . PDN=SENDWKS . I O . OWN* . 

* , 

REWIND. DN=XYZ. 

REWIND, DN=Q. 

• . 

GETP3D . OXYZ.OXOUT1 . DATATYPE—XYZ . DSSK IP—100 . NUMDS*1 50 . 

GETP3D , I=Q,O=Q0UT 1 . DATATYPES . DSSK I P=1 00 . NUMDS=1 50 . 

• . 

SENDWKS , Dte=XOUT 1 , VDN= ' FML I RIS1 " j i a a password" : : **/u/ j i aa/pagan/38x .iris"', NOR EC 
SENDWKS , DN=OOUT 1 , VDf^= ’ FML I RI SI " j i a a pas swa rd" : : "/u/ j i aa/pagan/38q .iris"', NOR EC 

• . 

•. FETCH, DN=J0B2. TEXT* • [PAGAN]SEND. JCL: 2 ' . 

•.REWIND, DN=J0B2 . 

• .SUBMIT, DN=J0B2. 
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4. SINE.TXT 



JOB, JN-FLYNAVY.T-30. ERIC PAGENKOPF , X4269 . 
ACCOUNT , AC- , US- . UPW- . 

• . 

ACCESS , DN-SENDWKS . PDN^SENDWKS .10—. OWN- . 



CFT.ON-A.OFF-S. 

• . 

ASSIGN. DN-X, A-FT50. 

ASSIGN. DN-O. A=FT 60 . 

• . 

LDR . MAP-PART , SET-ZERO . 

• . 

SENDWKS . DN=0 . VDN= ' FML IR I SI " j ioo PASSWORD": :/u/j i oo/pogon/qs i n ins' .NOREC. 
SENDWKS, DN— X . VDN= ’ FML I R I SI " j ioo PASSWORD": :/u/j i oo/pogon/snk r05s ’ .NOREC, 
CSF . SAVE . CSFACCT— . CSFPSWD- . DN=X . PDN-SNKR05S . I D-FLYNAVY . 

• . 

/EOF 

PROGRAM POINTER 



THIS PROGRAM GENERATES PL0T3D FILES WHICH WILL PROOUCE AN 
AOA POINTER WHEN FUNCTION 0 IS SELECTED. XYZ FILES ARE SAVED 
ON CSF FOR FUTURE ACCESS VIA PROGRAM GETSIN.JCL. A SINGLE Q 
FILE OF PROPER DIMENSION IS SENT DIRECTLY TO THE IRIS WITH 
FILENAME QSIN. IRIS. 

VARIABLES: 

NPLOTS: TOTAL NUMBER OF PLOTS IN ONE CYCLE. 

PNLOCI: INITIAL INTEGER LOCATION OF POINTER 
PNLOC : INTEGER LOCATION OF POINTER. 

XSCALE: AXIS SCALING FACTOR 

YSCALE : SINE WAVE SCALING FACTOR 

XPOSIT: X-POSITIONAL FACTOR (SCREEN LOCATION) 

YPOSIT: Y— POSITIONAL FACTOR (SCREEN LOCATION) 



C 

INTEGER PNLOC. GD I M, PNLOC I 

DIMENSION SI NX(3 , 500) , S INY( 3 , 500) , 01 ( 3 . 500) .02 ( 3 . 500) 
DIMENSION 03(3.500) ,04(3.500) 

DATA Q1 .02.03.04/6000*1 ./ 

••• INITIALIZE 
READ (5.1) 

READ (5.2) NPLOTS, PNLOCI .XSCALE. YSCALE. XPOSIT. YPOSIT 
PI=ATAN( 1 . )* 4 . 

GDIM = 3 
PNLOC - PNLOCI 

DO 10 L-1 .NPLOTS 

••* GENERATE PL0T3D POINTER FILE 

DO 20 N— 1 , NPLOTS+1 

XLOC - N*2 . *PI/NPLOTS 
SINY(I.N) - (0.0) *YSCALE + YPOSIT 
SINX(I.N) - (XLOC)*XSCALE + XPOSIT 
SINY(2,N) - (SIN(XLOC) )* YSCALE + YPOSIT 
SINX(2.N) - (XLOC) • XSCALE + XPOSIT 
20 CONTINUE 

PLOC = PNLOC* 2. *P I /NPLOTS 
S I NY ( 3 . 1 ) - (0.0) *YSCALE 4- YPOSIT 
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non 



SINX(J.I) - (PLOC) .XSCALE + XPOSIT 
SINY(3 , 2) - ( S I N(PLOC) ) .YSCALE + YPOSIT 
SI NX (3, 2) - (PLOC).XSCALE + XPOSIT 

WRITE(50) GDIM, NPLOTS 

WRITE(50)((SINX(I ,J) , 1-1 , 3) , J-1 .NPLOTS). 

((SI NY (I, J), 1-1 .3). J-1 .NPLOTS) 

PNLOC - PNLOC + 1 
I F (PNLOC . EO . NPLOTS) THEN 

PNLOC - PNLOC + 1 - NPLOTS 

END IF 
10 CONTINUE 

... GENERATE DUK44Y 0 FILE 

WRITE(60) GDIM. NPLOTS 

WRITE(60) GDIM. GDIM. GDIM. GDIM 

WRI TE(60) ((QI(I.J). 1=1 .3) .J-1 .NPLOTS), 

+ ( (Q2( I . J ) , 1-1 .3) . J-1 .NPLOTS). 

+ ( ( Q3 ( I , J ) . 1=1 .3). J-1 .NPLOTS), 

+ ((04(1, J). 1=1 .3) .J-1 .NPLOTS) 

C 

1 FORMAT (IX) 

2 FORMAT (1X.2I8.4F8.4) 

END 

/EOF 

NPLOTS: PNLOC I : XSCALE: YSCALE: XPOSIT: YPOSIT: 

294 222 .08 .1 -1.0 -.5 



C 

+ 

C 
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CPPLOT.TXT 



JOB. JN-SKYHAWK.T-30. ERIC PAGENKOPF . X4269 . 

ACCOUNT . AC- . US- . UPW- . 

• . 

CFT.O^-A.OFF-S. 

• . 

CSF . ACCESS . CSFACCT- . CSFPSWD- . DN-GETP3D . PDN-GETP3D . I O-RFR I CC . 
ACCESS . DN-SENDWKS . PDN-SENDWKS , I D— STTRDM . OWN-RFTRDM . 

CSF. ACCESS. CSFACCT-. CSFPSWD- . DN-CPX 1 . PDN-SNKR05P. 1 0— FLYNAVY . 

CSF. ACCESS. CSFACCT-, CSFPSWD-. DN=CFX1 . PDN— SNKR05F , ID-FLYNAVY. 

CSF. ACCESS. CSFACCT-. CSFPSWD-. DN=01 . PDN-SNKR0SO. ID-FLYNAVY. 

REWIND. DN—CPX1 . 

REWIND. DN-CFX1 . 

REWIND. DN-01 . 

GETP3D . I -CPX 1 , 0-CPX . DATATYPE-XYZ . DSSK I P-49 . NUMDS-2 . 

GETP3D , I — CFX 1 .O-CFX. DATATYPE-XYZ. DSSK IP-49, NUMDS-2. 

GETP3D.I-01. 0=0 , DATATYPE-Q. DSSK I P-49 . NUMDS-2 . 

REWIND. DN-CPX 
REWIND . DN-CFX . 

REWIND.DN-Q. 

ASS I GN . DN-CPX . A-FT50 . 

ASS I GN . DN-CFX . A-FT60 . 

ASSIGN. DN=0. A— FT70 . 

• . 

LDR . MAP-PART . SET-ZERO . 



/EOF 

PROGRAM CPPLOT 



THIS PROGRAM READS STORED. TRUE VALUES OF XOC. CP AND CF. 
AND PRINTS THEM OUT TO TAPE 06. 

VARIABLES: 

NUMDS: NUMBER OF DATA SETS IN THE READ FILE 
NPLOT : PLOT NO. OF FIRST DATA SET (DSSKIP + 1) 

TAPE ASSIGNMENTS: 

TAPE 05: INPUT DATA 
TAPE 06: OUTPUT DATA 
TAPE 50: TRUE CP VALUES 
TAPE 60: TRUE CF VALUES 
TAPE 70: Q FILES 



DIMENSION CPX(3.49) ,CFX(3 . 49) ,CPY(3 . 49) ,CFY(3 . 49) 
DIMENSION Q1 (161 .41 ) .02(161 .41) .03(161 .41 ) ,04( 161 .41 ) 
DIMENSION K0DE(4) . LINE(90) 

DATA K00E/1H , 1H+ , 1 HI , 1 H»/ 

... READ INPUT DATA k INITIALIZE 

READ (5.1) 

READ (5.2) NLAOS . NPLOT 
DO 9 1=1 .90 

LINE(I) - KODE( 1 ) 

CONTINUE 

... READ k SAVE TRUE VALUES 
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CM K> * O 



1-1 , IMAX) . J-1 , KMAX ) . 
1-1 . IMAX) .J-1 .KMAX) 



1 1 



10 

C 

1 



DO 10 LI-1 .NUMDS 

READ (50) IMAX. KMAX 
READ ( 50) ((CPX(I.J). 

+ ((CPY(I.J). 

READ( 60) IMAX. KMAX 
READ( 60) ((CFX(I.J), 1-1 . IMAX) . J-1 .KMAX) . 

+ ((CFY(I.J), 1-1 . IMAX) .J-1 .KMAX) 

READ(70) IMAXQ . KMAXO 
READ( 70) AMINE.ALFAD.REYREF.TIME 
READ(70) C FQ1 ( I . J). 1-1 . IMAXQ) . J-1 , KMAXQ) . 
+ ( v 02 C I . J ) . 1-1 . IMAXQ) .J-1 .KMAXQ) , 

+ ((03(1. J). 1-1 . IMAXO) .J-1 .KMAXQ) . 

+ ( (Q4( I . J ) , 1-1 , IMAXQ) .J-1 .KMAXQ) 

+ .2 • AMINE*. 2) »• 3.5 - 1 . 

7 • AMINF..2) 

+ 4.5 



(CP0 - CPY ( 3 . L) ) + 4.5 
(CP0 —CPY ( 2 . L) ) + 4.5 



CP0 - (1 
CP0 - CP0 / ( . 

K0 = 30. * CP0 
DO 11 L-1 .KMAX 
K — 30. • 

K 1 - 30. « 

I F(K . LT . 1 ) K - 1 
I F ( K 1 . LT. 1) K 1 - 1 
I F(K . GT . 90) K - 90 
I F ( K 1 .GT. 90) K1 - 90 
LINE(K0) - KODE(3) 

LINE(K) = KODE( 2) 

LINE(K1 ) - K0DE(4) 

I F( L . EQ . 1 )THEN 

WR I T E ( 6 . • ) ' PLOT AOA TIME MACH REY NO. * 

WRITE(6.3) NPLOT . ALFAD , T IME , AMINE .REYREF 
WRITE(6, *) ’ XOC CFL CFU CPL CPU * 

WRITE(6.4) CPX(1 ,L).CFY(2.L),CFY(3,L) ,CPY(2, L) ,CPY(3, L) . LINE 
ELSE 

WRITE(6.4) CPX(1 ,L).CFY(2.L) . CFY(3 . L) . CPY(2 . L) . CPY(3 . L) . LINE 
END IF 

LINE(K1 )-KODE( 1 ) 

LINE(K) -KODE( 1 ) 

CONTINUE 

NPLOT - NPLOT + 1 
WRITE(6. 1 ) 

CONTINUE 

FORMAT (IX) 

FORMAT (IX. 217) 

FORMAT (IX. I4.3F9.4.F10.0) 

FORMAT (IX. F6.4.4F8.4.90A1 ) 



END 

/EOF 

NUMDS: NPLOT: 
2 50 
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6. GETX.F 



program getx 

c 



THIS PROGRAM READS COMBINED (MULTIPLE DATA SET) FILES 
OF XYZ AND 0 PLOTTING DATA AND SEPARATES THEM INTO 
INDIVIDUAL FILES FOR PLOT3D SUBMITTAL. 

VARIABLES: 

CFXNAME : NAME OF COMBINED XYZ FILE 
CFONAME: NAME OF COMBINED 0 FILE 

FXNAME: NAME Of INDIVIDUAL X FILES. CHAR VARIABLE 
FONAME: NAME OF INDIVIDUAL 0 FILES. CHAR VARIABLE 
DSSKIP: NUMBER OF DATA SETS TO SKIP WHEN 
READING COMBINED FILE 
DSREAD : NUMBER OF DATA SETS TO READ 
INTVL: INTERVAL BETWEEN DATA SETS SENT TO FILE 



C 

INTEGER DSSKIP. DSREAD. COUNT 

CHARACTER FXNAME* 1 2 . FONAME* 1 2 . CFXNAME* 10 , CFONAME* 1 0 
DIMENSION X(250, 100) ,Z(250. 100) 

DIMENSION 01(250.100) .02(250. 100) ,03(250, 100) .04(250. 100) 
C 

••• USER INPUT: 

C 

DATA XSCALE.YSCALE/1 . . 1 ./ 

DATA XPOS I T . YPOS I T /— 1 . , - 1 . / 

DATA DSSKIP/10/.DSREAD/3/. INTVL/1/ 

CFXNAME - ’ 38x . iris’ 

CFONAME - ' 38q . iris’ 

C 



c 



c 



c 



10 

c 



OPEN (UNIT-91.FI LE=CFXNAME . STATUS— ’ OLD ’ . FORM- ’ B I NARY ’ ) 
OPEN(UNI T-92 , FI LE-C FONAME . STATUS- ’OLD ’ .FORM- ’BINARY’ ) 

FXNAME « ’ x000. iris’ 

FONAME = ’ q000 . iris’ 

COUNT - 0 



IF (DSSKIP. GT.0) THEN 

DO 10 ISKIP = 1 .DSSKIP 



READ(91 ) 
READ(91 ) 

READ(92) 

READ(92) 

READ(92) 



END IF 



CONTINUE 



I MAX , KMAX 

( (X( I ,K) . 1-1 . IMAX) ,K«1 .KMAX) , 
((Z(I.K).I-I . IMAX ) ,K=1 .KMAX) 
IMAX. KMAX 
A.B.C.D 

( (01 (I ,K) . 1=1 , IMAX) .K-1 .KMAX) . 
( (02(1 ,K) . 1-1 . IMAX) , K— 1 .KMAX) , 
((03(1 .K). 1-1 . IMAX) .K-1 .KMAX) , 
((04(1 ,K) . 1=1 . IMAX) .K-1 .KMAX) 



DO 20 IREAD - 1 .DSREAD 
COUNT « COUNT + 1 
WRITE ( FXNAME ( 2 : 4) . 1 00) IREAD 
WRITE ( FONAME ( 2 : 4) , 1 00) IREAD 
READ(91) IMAX. KMAX 

READ(91 ) ( (X( I .K) . 1 — 1 . I MAX). K-1 .KMAX) . 
+ ((Z(I ,K). 1-1 . IMAX), K-1 .KMAX) 
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12 
1 1 



20 

C 



100 



IF (COUNT . EO. INTVL) THEN 

OPEN ( UN I T« 1 . F I LE-FXNAME . FORM- 1 B I NARY ’ ) 
DO 1 1 1-1 . IMAX 

DO 12 K-1 , KMAX 

X( I ,K)-X( I , K) -XSCALE+XPOS IT 
Z( I ,K)-Z( I , K)*YSCALE+YPOSI T 
CONTINUE 
CONTINUE 

WRITE(I) IMAX, KMAX 

WRITE( 1 ) ( (X( I ,K) , 1-1 , IMAX) ,K=1 , KMAX ) , 
+ ( (Z( I ,K) , 1-1 , IMAX) . K-1 , KMAX) 

CLOSE(I) 

END IF 

READ(92) IMAX, KMAX 
READ( 92) A.B.C.D 

READ( 92) ( (01 ( I , K) . 1 = 1 . I MAX). K-1 .KMAX) , 

+ ((Q2( I ,K) , 1-1 , IMAX) ,K=1 .KMAX) , 

+ ( (Q3( I ,K) , 1-1 , IMAX) ,K=1 .KMAX) . 

+ ((Q4( I ,K) , 1-1 , IMAX) ,K-1 .KMAX) 

IF (COUNT. EQ. INTVL) THEN 

OPEN ( UN I T-2 . F I LE-FONAME . FORM- ’ B I NARY 1 ) 
WRITE(2) IMAX, KMAX 
WRITE(2) A.B.C.D 

WRITE(2) ( (Q1 ( I ,K) . 1=1 . IMAX) , K-1 .KMAX) , 
+ ((02(1 ,K) . 1-1 . IMAX) , K-1 .KMAX) , 

+ ((03( I ,K) , 1-1 . IMAX) ,K-1 .KMAX) . 

+ ((04( I ,K) . 1-1 . IMAX) .K-1 .KMAX) 

CL0SE(2) 

COUNT - 0 

END IF 
CONTINUE 

CL0SE(91 ) 

CLOSE(92) 

FORMAT( 13.3) 

STOP 

END 
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7. GETCP.F 



PROGRAM GETCP 
C 



THIS PROGRAM ACCESSES COM3INED CP OR CF FILES ON THE 
IRIS ACCOUNT, SEPARATES THEM INTO INDIVIDUAL PLOTTING 
FILES AND SCALES THEM TO MEET SPECIFIED REQUIREMENTS. 

A DUMMY 0 FI L r IS ALSO GENERATED (Q.IRIS). 

VARIABLES: 

CPNAME : COMBINED FILE NAME 
FNAME: OUTPUT FILE NAME, CHARACTER VARIABLE 
DSSKIP: NUMBER OF DATASETS IN THE COMBINED FILE 
TO SKIP 

DSFILE: NUMBER OF DATASETS IN THE COMBINED FILE 
TO SEPARATE AND FILE 
XSCALE : AXIS SCALING FACTOR 

YSCALE: CP/CF SCALING FACTOR (RELATIVE MAGNITUDE) 
XPOSIT: X POSITION FACTOR (SCREEN LOCATION) 
YPOSIT: Y POSITION FACTOR (SCREEN LOCATION) 



C 

INTEGER DSSKIP. DSFILE 
CHARACTER FNAME* 1 1 , QNAME*6 , C FNAME* 10 
DIMENSION X(3.75) ,Z( 3.75) ,XS( 3,75). ZS(3. 75) 
DIMENSION 01(3,75) .02(3,75) .03(3,75) ,04(3,75) 
DATA 01,02,03.04/900*1./ 

C 

• •• USER INPUTS: 

C 

DATA DSSKIP/250/.DSFI LE/44/ 

DATA XSCALE/. 5/. YSCALE/-. 05/ 

DATA XPOS I T/- 1 . 0/ . YPOS I T/ . 5/ 

CFNAME - 'snkr05p’ 

FNAME - 'epee©. iris' 

C 



C 

OPEN ( UN I T-90 , F I LE=CFNAME . STATUS- ’ OLD ’ . FORM- ’ 81 NARY ’ ) 
C 

QNAME - ' q . i r i s ’ 

C 

••• SKIP DSSKIP FILES 
C 

IF (DSSKIP. GT 0) THEN 

DO 10 ISKIP - 1 .DSSKIP 
READ(90) I MAX . KMAX 

READ(90) ((X(I ,K) . 1-1 , IMAX) ,K«1 .KMAX) . 

+ ((2(1 ,K). 1-1, IMAX). K-1, KMAX) 

10 CONTINUE 

END IF 
C 

••• SEPARATE AND SCALE DSFILE FILES 
C 

DO 40 IFILE - 1 .DSFILE 

WRI TE( FNAME(3 : 5) .100) IFILE 

OPEN ( UN I T- 1 . F I L E-FNAME . FORM* ' B I NARY ' ) 

READ(90) IMAX. KMAX 

READ(90) ((X( I ,K) . 1-1 . IMAX) .K-1 .KMAX) . 

+ ((Z(I ,K) , 1-1 . IMAX) .K-1 .KMAX) 
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DO 30 I - 1 . I MAX 

DO 20 K - 1 , KMAX 

XS( I ,K)-X( I ,K)»XSCALE + XPOSIT 
ZS ( I ,K)=Z(I ,K)»YSCALE + YPOSIT 
20 CONTINUE 

30 CONTINUE 

WRITE(l) I MAX .KMAX 

WRITE(1 ) ((XS( I ,K) , 1*1 . I MAX) , K*1 .KMAX) , 

+ ( (ZS ( I ,K) , 1-1 . IMAX) ,K*1 , KMAX) 

C LOS E ( 1 ) 

40 CONTINUE 
C 

GENERATE DUfcMY 0 FILE 
C 

OP EN ( UN I T«= 1 . FI LE«=ONAME . FORM= ’ B I NARY 1 ) 

WRITE(I) IMAX , KMAX 

WRITE(I) IMAX, IMAX, IMAX, IMAX 

WRI TE( 1 ) ((01 ( I ,K) , 1-1 . IMAX) ,K=1 .KMAX) . 

+ ((Q2(I .K) . I«1 , IMAX). K»1 .KMAX), 

+ ( (03( I ,K) . 1=1 , IMAX) ,K-1 .KMAX) . 

+ ((Q4( I ,K) . 1*1 . IMAX) ,K*1 .KMAX) 

CLOSE( 1 ) 

C 

100 FORMAT (13.3) 

C 

STOP 

END 
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8. GETS1N.F 



PROGRAM GETSIN 
C 



THIS PROGRAM READS A COMBINED FILE OF AOA POINTER 
DATA (GENERATED WITH PLOT3D FUNCTION 0) AND SEPARATES 
IT INTO INDIVIDUAL PLOTTING FILES. 

VARIABLES: 

CFNAME : NAME OF COMBINED FILE 

FNAME : NAME OF INDIVIDUAL FILES. CHAR VARIABLE 
DSSKIP: NUMBER OF DATA SETS TO SKIP WHEN 
READING COMBINED FILE 

DSFILE: NUMBER OF DATA SETS TO SEND TO INDIVIDUAL 
FILES 



C 

INTEGER DSSKIP. DSFILE 
CHARACTER FNAME* 12. CFNAME* 10 
DIMENSION X(3, 294) .2(3.294) 

C 

••• USER INPUT: 

C 

DATA DSSKIP/250/.DSFI LE/44/ 
CFNAME = ’snkr05s’ 

C 



C 

OPEN(UNIT=90.FILE=CFNAME.STATUS=’OLD' , FORt^'BINARY’ ) 
C 

FNAME = ’ s i n000 . i r i s ’ 

C 

IF (DSSKIP. GT.0) THEN 

DO 10 I SKIP - 1 .DSSKIP 
READ( 90) I MAX , KMAX 

READ( 90) ((X( I . J) . 1 = 1 . I MAX) ,J=1 .KMAX) . 

+ ( (Z ( I .J) , 1=1 . IMAX).J=1 .KMAX) 

10 CONTINUE 

END IF 
C 

DO 20 IFILE = 1 .DSFILE 

WRITE (FNAME(4:6) , 100) IFILE 

OPEN( UN I T= 1 , F I LE= FNAME , FORV^ * B I NARY * ) 

READ(90) I MAX . KMAX 

READ(90) ((X(I .J) , 1=1 . I MAX) , J=1 .KMAX). 

+ ((Z(I,J). 1=1. IMAX),J=1. KMAX) 

WRITE(I) I MAX, KMAX 

WRITE(I) ((X(I.J). 1=1 . I MAX ) , J=1 .KMAX), 

+ ( (Z( I , J) . 1 = 1 . IMAX) ,J=1 .KMAX) 

CLOSE(I) 

20 CONTINUE 
C 

CLOSE(90) 

100 FORMAT( 13.3) 

STOP 

END 
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C. PROCEDURAL DOCUMENTATION 



After proper initialization (user-inputs) as described above, all codes or JCL s \vli'u.h 
produce or transfer Plot3D files aie submitted to the Cray utilizing standard CSUB 
commands. Once the data is resident on the IRIS and the user is logged on and in the 
proper directory level, the following procedures may be utilized for graphics generation. 

1. Separate the combined file into discrete data sets. 

Edit "user inputs" section of GE i’X.E (or other files as appropriate) to define scal- 
ing (ultimate location on the display), and identify the combined file to access and 
those data sets to be separated: 

vi getx.f 

Compile the editted program: 

ill getx.f 

Run the compiled version by entering the filename (a. out by default). 

a. out 

Discrete data sets will appear on the account with names corresponding to those 
in the PlotiD initialization files. Check by using the Is command. 

2. Use Plot3D to generate .gra files. 

Use the mv commands to give the proper initialization file tire filename 
"plot3dini.com". (The IRIS will not save multiple versions of files with the same 
filename. Instead, older versions of the file will be deleted. In order to avoid in- 
advertent deletion of initialization files, multiple mv commands may be required.) 

inv xini.com plot3dini.com 

If processing How field files (X and Q), edit the initialization Hie to identify desired 
function and number of contours. This involves changing only two lines at the top 
of the file. 



vi plot3dini.com 

Create .gra files. 

plot3x 

Rename the resultant combined graphics file. 

mv qOOl.gra somefilnm.gra 

After creation of as many .gra files as required from the separated files, dean up 
the account using the following wildcards. 

rm q " .iris 

rm x*.iris 
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. Interactive viewing. 

Run the Graphics Animation System software, 
gas 

Input .gra files to GAS by selecting on the main menu: 

ARCGRAPM lile input 

On the submenu, select: 

load entire file 

In response to prompts, enter the .sequential object number and file name 
("someflnm.gra"). For color plots, press RETURN when prompted for the 
colormap. 

View the data by selecting on the main menu: 

view data 

Utilize the mouse to manipulate the display. 
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APPENDIX B. SANKAR NAVIE STOKES SOLVER 



(JCL) 

• . 

/EOF 

PROGRAM MAIN 



SANKAR NAVIER STOKES CODE: SOLVES TWO DIMENSIONAL FLOWS PAST 

ARBITRARY GEOMETRIES USING ADI PROCEDURE. THIS VERSION SAVES 
MULTIPLE PLOT3D DATA SETS IN SINGLE FILES. 

VARIABLES: 

TITLE: (60 CHARACTERS MAX. ) 

IMAX: NUMBER OF X COORD LOCATIONS 
KMAX : NUMBER OF Y COORD LOCATIONS 
ITEL: GRID POINT AT LOWER TRAILING EDGE 
I TEU : GRID POINT AT UPPER TRAILING EDGE 
DT: SIZE OF TIME STEP 

WW: EXPLICIT ARTIFICIAL VISCOSITY TERM 
ALFA: MEAN ANGLE OF ATTACK (DEGREES) 

ALFA1 : AMPLITUDE OF OSCILLATION (DEGREES) 

ALFAI: ANGLE BELOW WHICH MODIFIED TURBULENCE MODEL USED TO 
COMPUTE EDDY VISCOSITY 
REDFRE: REDUCED FREQUENCY 
AM INF: FREE STREAM UACH NUMBER 

REYREF: REYNOLDS NUkCER (MILLIONS: NEG. = INVISCID FLOW) 

DNMIN: DISTANCE OF FIRST POINT OFF OF WALL 

TSTART : TIME FOR CALCULATIONS TO START ON PRESENT RUN 

FORMAT: OLDSLN FORMAT ; 3 . 0-P LOT3D FI LES . -3 . 0=0 MATRICES 

RSTRT: TRUE = STORED DATA READ TO CONTINUE ITERATION 

PITCH: TRUE - AIRFOL AOA OSCILLATES 

PLUNGE: TRUE = AIRFOIL OSCILLATES VERTICALLY 

FNU: NUMBER OF UPPER SURFACE DEFINITION POINTS 

FNL : NUMBER OF LOWER SURFACE DEFINITION POINTS 

FSYM : SYWETRY FLAG (1 = SYMMETRIC) 

X: AIRFOIL GEOMETRY DEFINITION POINTS 
Y: AIRFOIL GEOMETRY DEFINITION POINTS 
CSTP : NUMBER OF STEPS COMPLETED 
CPLT : NUMBER OF DYNAMIC PLOTS COMPLETED 

NSTP : NUMBER OF TIME STEPS TO BE COMPLETED ON PRESENT RUN 
PSTP : NUMBER OF TIME STEPS BETWEEN PLOT DATA OUTPUT 

TAPE DEFINITIONS: 

TAPE 05: INPUT DATA 
TAPE 06: OUTPUT DATA 

TAPE 07: FLOW FIELD INPUT DATA FOR RESTARTS 
TAPE 20: PL0T3D XYZ FILE STORAGE 
TAPE 21: PL0T3D 0 FILE STORAGE 
TAPE 50: PLOT 3D CP XYZ FILE STORAGE 
TAPE 60: PL0T3D CF XYZ FILE STORAGE 
TAPE 90: PREVIOUS RUN X FILE STORAGE 
TAPE 91: PREVIOUS RUN 0 FILE STORAGE 
TAPE 92: PREVIOUS RUN CP FILE STORAGE 
TAPE 93: PREVIOUS RUN CF FILE STORAGE 



INTEGER PSTP.CSTP.CPLT 
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c 

c • • • 

c 



c 

c • • • 

c 



c 

c 

c 



c 

c 

c 



7 

c 

c • • • 

c 



COWON/SUR F/PSUR (161) 

COMADN/ FIX/OMEGA 
CO*ADN/MUTUR/CMU( 161 ,41 ) 

CO»ADN/SK I NC F/C F ( 1 6 1 ) 

CC***DN/GRID1/X(161 ,41). 2(161,41) 

CO*ADN/PAR/GA>A4A . REYREF . ALFA , ALF A 1 , REDFRE . AM I NF . ALFA I 
CO440N/0GR I D/DT . I MAX , KMAX , I TEL, ITEU 
C0U40N/GR I O/YACOB (161,41) 

COMMON /DAMP/WW , WW I 

CCA4ADN/FL0W/Q1 ( 161 .41 ) .02(161 . 41 ) ,Q3( 161 , 41 ) ,Q4( 161 ,41 ) 

DIMENSION TITLE( 15) ,TYTLE(15) ,ALPHA(96) ,CPTH(97 . 96) . XTH(97) 
DIMENSION XR( 161 .49) ,ZR(161 .49) ,01R(161 .41) ,02R(161 ,41) 

DIMENSION 03R (161 .41 ) ,Q4R(161 .41 ) 

COM*DN/LOC 1 C/RSTRT .PITCH .PLUNGE 
LOGICAL RSTRT .PITCH. PLUNGE 

READ INPUT DATA 

READ (5.5) 

READ (5,1) TITLE 
READ (5.5) 

READ (5.2) I MAX. KMAX. DT.WW. ALFA. ALFA1 .ALFAI .REDFRE. AMINF 
READ (5,5) 

READ (5.3) ITEL. ITEU. REYREF. DNMI N . TSTART . FORMAT. RSTRT. PITCH, PLUNGE 
READ (5.5) 

READ (5.4) CSTP . CPLT . NSTP . PS TP 

INITIALIZE 

GAMMA =1.4 
PI = ATAN( 1 . )*4. 

SET ALFAD FOR STEADY STATE PL0T3D OUTPUT 
ALFAD = ALFA 

FORCE DT TO BE EQUAL TO UNITY FOR STEADY FLOW PROBLEMS 

THIS INVOKES THE RELAXATION MODE 

IF(REDFRE. LE. 0.001 ) DT = 1.0 

REYREF = REYREF * 1 . E+06 

NITN = CSTP + NSTP 

NPLOTS = CPLT 

CPLT - CPLT + 1 

CSTP = CSTP + 1 

DENSITY NORMALISED WITH RESPECT TO ROINF 

VELOCITIES NORMALISED WITH RESPECT TO AlNF 

TOTAL ENERGY NORMALISED WITH RESPECT TO (ROINF. AI NF* A INF) 

TOTEN=AMI NF.AMI NF *0.5+1 . / ( GAK44A* (GAJJ4A— 1 . )) 

ALFA = ALFA • ATAN( 1 . ) / 45. 

ALMEAN = ALFA 

ALFAI « ALFAI • ATAN( 1 . ) / 45. 

ALFAI = ALFAI • ATAN( 1 . ) / 45. 

ALFAS = ALMEAN - ALFAI 

WR I T E ( 6 ) ' PLOT I TN TIME AOA CL CD CM’ 

UINF = AMINF 
VINF = 0.0 
DO 7 1=1 . IMAX 
DO 7 K=1 .KMAX 
01(1 ,K)=1 . 

02 ( I ,K)=UINF 
Q3( I , K )=V I NF 
Q4( I , K)=TOTEN 
CONTINUE 

INPUT STORED FLOW FIELD DATA 

IF(RSTRT)THEN 

I F (FORMAT. LT. 0.0) THEN 

READ(7) TIME. 01 .02.03.04 

ELSE 

READ(7) IMAX. KMAX 
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oooo - ooo 



+ 

+ 



READ( 7 ) AMINF , ALFAD . REYREF .TIME 
READ( 7 ) ((01(1, J), 1=1 , I MAX) , J=1 , KMAX ) , 

((02(1, J). 1=1 , I MAX) , J=1 , KMAX) , 
((03(1, J), 1=1 .IMAX). J=1 .KMAX), 
((04(1, J), 1=1 , IMAX) . J— 1 .KMAX) 



END IF 
IF(NPLOTS.GT.0)THEN 
DO 9 K 1=1 , NPLOTS 



READ (90) 
READ(90) 

WRITE (20) 
WRITE (20) 



CONTINUE 
DO 10 K2=1 .NPLOTS 



IMAXR.KMAXR 
( (XR( 1 . J) , 1 = 
((ZR(I.J). I- 
IMAXR.KMAXR 
((XR(I.J). I. 
( (ZR( I . J) , I- 



1 , IMAXR) , Ja 
1 . IMAXR) , J= 



1 . KMAXR) . 
1 .KMAXR) 



1 . IMAXR) . J=1 .KMAXR) , 
1 . IMAXR) . J=1 .KMAXR) 



+ 

+ 

+ 



+ 

+ 

+ 



READ (91 ) 
READ(91) 
READ(91) 



WRITE(21 ) 
WR I TE(21 ) 
WRITE(21 ) 



IMAXR.KMAXR 

AMINFR . ALFADR .REYREFR , TIMER 
((OIR(I.J). 1=1 . IMAXR) ,J=1 .KMAXR) 
1=1 . IMAXR) . J=1 , KMAXR) 
1=1 . IMAXR) . J=1 .KMAXR) 
1=1 , IMAXR) , J=1 .KMAXR) 
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( (Q2R( I . J) . 

( ( 03R ( I . J ) . 

( (Q*R( I . J ) . 

IMAXR.KMAXR 
AMINFR. ALFADR 
((QiR(l.J). 

( (Q2R( I , J) , 

( ( 03R ( I .J), 

( (Q4R( I . J) , 



REYREFR, TIMER 
1 = 1 . IMAXR) , J=1 .KMAXR) 
1 = 1 . IMAXR) . J«=1 .KMAXR) 
1=1 . IMAXR) . J=1 .KMAXR) 
1=1 . IMAXR) . J=1 , KMAXR) 



1 1 
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CONTINUE 
DO 11 K3=1 , NPLOTS 

READ(92) IMAXR.KMAXR 

READ(92) ((XR(I.J), 1=1 .IMAXR) ,J=1 .KMAXR) , 
((ZR(I.J). 1=1 . IMAXR) ,J=1 .KMAXR) 
WRITE(50) IMAXR.KMAXR 

WRITE(50) ((XR(I.J), 1=1 .IMAXR), J=1 .KMAXR), 
( (ZR( I . J) . 1 = 1 , IMAXR) ,J=1 .KMAXR) 

CONTINUE 

DO 12 K4=1. NPLOTS 

READ(93) IMAXR.KMAXR 

READ(93) ((XR(I.J), 1 = 1 . IMAXR) . J = 1 .KMAXR) . 

((ZR(I.J), 1=1 .IMAXR) . J= 1 .KMAXR) 
WRITE (60) IMAXR.KMAXR 

WRITE (60) ((XR(I.J), 1=1 .IMAXR) ,J=1 .KMAXR) . 

((ZR(I.J), 1=1 . IMAXR) ,J=1 .KMAXR) 

CONTINUE 



END IF 

END IF 

I F( TSTART . GE . 0 . ) 
IF( .NOT. (RSTRT)) 



TIME = TSTART 
TIME = 0. 



••• GENERATE COMPUTATIONAL GRID, DEFINE SURFACE COORD O 0 AOA , 
ROTATE TO INITIAL AOA AND COMPUTE METRICS 

CALL AIRFOL( IMAX . KMAX . I TEL , I TEU) 

IF (REYREF. GT.0) CALL CLUSTR(DNMIN) 

TEDGE = X ( I TEL + 48.1) 

DO 15 I = 1.97 

XTH(I) = X ( I+ITEL-1 . 1 )-TEDGE 
5 CONTINUE 

CALL ROTGRID(X,Z. IMAX. KMAX. ALFAS) 

CALL METRIC 

... ITERATIVE LOOP 



DO 1000 I TN=1 , NSTP 
TIME = TIME + DT 
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<>•• ROTATE GRID TO NEW ANGLE OR SET TO CORRECT ANGLE FOR RESTARTS 
C 

IF (PITCH) THEN 

OMEGA - 2. • REDFRE .AMINF. SIN(REDFRE • 2. • TIME • AMINF) 
1 *ALFA1 

ALOLD - ALMEAN - ALFA1 • COS(2. • REDFRE . AMINF • 

1 (TIME - DT)) 

ALFA - ALMEAN - ALFA1 • COS(REDFRE • 2. • TIME » AMINF) 
ALFAD = ALFA • 45. / ATAN( 1 . ) 

DALFA = ALFA - ALOLD 

I F ( RSTRT . AND . I TN . EO . 1 ) DALFA = ALFA - 2.*ALFAS 
CALL ROTGRID(X.Z, IMAX. KMAX, DALFA) 

CALL METRIC 
END IF 

IF (PLUNGE) THEN 

OMEGA = 2. • REDFRE • AMINF 
ALFA - ALMEAN 
END IF 

C 

C* • • COMPUTE THE SOLUTION BY ADI TECHNIQUE 
C 

CALL SLPS( ITN. OMEGA. DALFA) 

C 

C... APPLY BOUNDARY CONDITIONS 
C 



C 

C* • • 

c 



c 

c • * * 

c 



1000 

c 

1 

2 

3 

4 

5 

6 
C 

C 

c»». 

C 



CALL WALLBC 

GENERATE PLOT 3D FILES 

IF(CSTP.EQ. (CPLT.PSTP)) THEN 
WRITE(20) I MAX , KMAX 

WRITE(20) ((X(I,J). 1=1,1 MAX ) , J=1 . KMAX ) , 

1 ((Z(I,J). 1=1,1 MAX ) , J=1 . KMAX ) 

WRITE(21) I MAX , KMAX 
WRITE(21 ) AMINF, ALFAD. REYREF, TIME 
WR I T E ( 2 1 ) ((QI(I.J), 1 = 1 , IMAX) , J=1.KMAX), 

1 ((Q2(I. J). 1=1, IMAX). J=1 , KMAX ) , 

2 ((03(1. J). 1=1, IMAX), J=1 , KMAX ) . 

3 ((04(1. J). 1=1, IMAX), J=1 , KMAX ) 

GENERATE PERFORMANCE COEFFICIENTS 

CALL LOAD(PSUR , CL , CD , CM , ALFAS) 

AOA = ALFA. 180. /PI 

WR I T E ( 6 , 6 ) CPLT , CSTP , T I ME , AOA ,CL ,CD , CM 
CALL CPPLOT ( I TEL , I TEU , AMINF . X( 1 , 1 ) . Z( 1 . 1) ,PSUR) 
CPLT = CPLT + 1 
END IF 

CSTP = CSTP + 1 
CONTINUE 

FORMAT (1X.15A4) 

FORMAT (1X.2I7.7F8.4) 

FORMAT ( 1X.2I7.4F8.5.3L7) 

FORMAT (IX, 4 I 7) 

FORMAT (IX) 

FORMAT (IX. I3.I6.F8.3.F9.5.3F8.4) 

END 



SUBROUTINE AMAT 1 ( K , IMX1.XIX.XIZ.XIT) 

CCA440N/ FLOW/01 (161 . 41 ) . Q2( 1 61 . 41 ) . Q3( 161 . 41 ) .Q4( 1 61 , 41 ) 
CCA440N/PERTR/DQ1 (161 .41 ) .DQ2(161 .41 ) .003(161 .41 ) ,D04( 161 .41 ) 
CCA440N /AM/ A ( 4 , 4 . 161) 

C0fcl40N/PAR/GAK*4A. REYREF. ALFA. ALFA 1 . REDFRE . AMI NF , ALFAI 
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c 

c • • • 

c 



DIMENSION XIT(16l.4l),XIX(161,41).XIZ(161.4l) 

REAL K0.K1 , K2 

AMAT1 COMPUTES THE COEFFICIENT MATRIX DE/DQ DURING XI SWEEP 



1000 



C 

C 

C 



GM1 

DO 1000 I 

K0 

K1 

K2 

U 

W 

EBYR 
PHI 2 
THETA 
A(1 ,1 .1) 
A(1 .2.1) 
A(1 .3.1) 

A ( 1 . 4 , I ) 

A ( 2 , 1 . I ) 
A( 2 , 2 . I) 
A(2 . 3 , I ) 
A(2 . 4 . I ) 
A(3 , 1.1) 
A(3 , 2 . I ) 
A(3 . 3 . I ) 
A(3 . 4 . I ) 
A( 4 , 1 , I ) 
A(4.2. I) 
A(4 . 3 , I ) 
A(4 , 4, I ) 
CONTINUE 
RETURN 
END 



- GAMMA - 1 . 

« 2 . IMX1 

- X I T ( I .«) 

- X I X( I ,K) 

- XI Z( I ,K) 

* 02(1. K) / 01(1, K) 

= 03(1 ,K) / 01(1 . K) 

= 04 ( I , K) / 01 (I .K) 

- 0.5 • GM1 • (U • U + W • W) 



* 


K1 


• 


U + K2 


• W 






- 


K0 












* 


K 1 












= 


K2 












m 


0 












m 


K 1 


• 


PH I 2 — 


U * 


THETA 


rr 


K0 


+ 


THETA - 


- K 1 


0 


(GM1 


« 


K2 


• 


U - GM1 


• K 1 


• W 


s 


K 1 


• 


GM 1 








c s 


K2 


• 


PH I 2 — 


W • 


THETA 




K 1 


• 


W - K2 


« GM 1 


• U 




K0 


+ 


THETA - 


K2 


• 


(GM1 




K2 


• 


GM1 









=■= THETA • (2. • PHI 2 - GAMMA • EBYR) 

- K1 • (GAMMA • EBYR - PHI2) — GM1 . U • THETA 

- K2 • (GAMMA • EBYR - PHI2) — GM1 • W • THETA 

« K0 + GAMMA • THETA 



SUBROUTINE AMAT2( I , KMX1 . ZETAX . ZETAZ . ZETAT) 

C0M40N/FL0W/Q1 (161 . 4 1 ) , 02 ( 1 61 , 4 1 ) . 03( 1 61 , 4 1 ) . 04( 1 61 . 41 ) 
COMAON/PERTR/DOI (161 ,41).D02(161 .41). 003(161 .41). 004(161 .41) 
COMAON/AM/A ( 4 . 4 . 1 6 1 ) 

COMMON/PAR/GAMMA , REYREF . ALFA , ALFA1 . REDFRE . AM I NF . ALFA I 
DIMENSION ZETAX (161 ,41) .ZETAZ (161 .41) ,ZETAT( 161 .41) 

REAL K0.K1 ,K2 



C 

C • • • 

c 



AMAT2 COMPUTES THE COEFFICIENT MATRIX DF/DQ DURING ETA SWEEP 



GM1 

DO 1000 K 

K0 

K 1 

K2 

U 

W 

EBYR 
PH 1 2 
THETA 
A ( 1 . 1 , K ) 

A ( 1 , 2 ,K) 
A( 1 , 3 , K ) 
A( 1 , 4 , K ) 
A(2 , 1 , K ) 
A(2 . 2 , K ) 
A(2 , 3 , K ) 
A(2 , 4 , K ) 
A(3 , 1 .K) 
A(3 , 2 . K ) 
A(3 . 3.K) 



GAMMA - 1 . 

2 , KMX1 
ZETAT(I ,K) 

ZETAX(I , K) 

ZETAZ(I .K) 

Q2( I ,K) / 01(1 . K) 

03(1 ,K) / 01 ( I , K ) 

04(1. K) / 01(1 .K) 

0.5 • GM1 • (U • U + W • W) 

K1 • U + K2 • W 

K0 

K1 

K2 

0 

K1 . PHI 2 -U • THETA 
K0 + THETA - K1 • (GM1-1.) • U 
K2 * U — GM1 • K1 • W 
K 1 • GM1 

K2 • PH 1 2 — W • THETA 
K 1 • W - K2 • GM1 • U 
K0 + THETA — K2 • (GM1-1 . ) • W 
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A(3,4,K) - K2 • GM1 

A( 4 , 1 , K) * THETA . (2. • PHI2 - GAM4A • EBYR) 

A(4.2,K) - K 1 . (GAM4A • EBYR - PHI2) - GM1 • U • THETA 

A(4.3,K) - K2 • (GAM4A • EBYR - PHI2) - GM1 • W • THETA 

A( 4 , 4 , K ) « K0 + GAM4A • THETA 

1000 CONTINUE 
RETURN 
END 
C 



C 

SUBROUTINE SLPS( I TN . OMEGA . DALFA) 

COMMON/ FLOW/01 ( 161 . 41 ) ,Q2( 161 .41 ) .03(161 .41 ) ,Q4( 161 .41 ) 

CCA440N/PERTR/DQ1 (161 ,41 ) ,DQ2( 161 ,41 ) .003(161 .41 ) ,004(161 ,41 ) 

COMMON /AM/ A (4, 4. 161 ) 

COMwtON/TR 10/00(4, 4, 161 .41). ^4(4, 4, 161, 41), EE(4, 4, 161, 41) 

1 . GG ( 4 , 161,41) 

COMMON/PAR/GAV44A . REYREF . ALFA . ALFA1 . REDFRE .AMI NF .ALFA I 
C0MM0N/DGR1D/DT. 1 MAX , KMAX . I TEL. ITEU 
COMMON/GR I O/YACOB (161,41) 

COMMON/DAMO/WW , WW 1 

COMMON/MTR I X/ X I X ( 1 61 . 41 ) . XI Z( 1 61 . 41 ) ,ZETAX(161 .41 ) , 

+ZETAZ( 161 .41 ) 

1 ,XIT(161 .41) ,ZETAT(161 .41) 

REAL 

DIMENSION DELTAT( 161 .41 ) 

C 

C* • • SUBROUTINE SLPS DOES THE BULK OF THE WORK FOR THE AC! ALGORITHM. 

C* • • IT CALLS FLUX AND COMPUTES RIGHT HAND SIDE DURING THE TWO SWEEPS. 

C* • • ASSEMBLES THE COEFFICIENT MARICES. ADDS IMPLICIT AND EXPLICIT 
C«*« DISSIPATION AND CALLS THE TRIDIAGONAL SOLVER TO OBTAIN THE FINAL 
C«»» SOLUTION. 

C ! ! ! ! ! SET VALUE OF IMPLICIT DAMPING COEFFICIENT 
WWI = 20.0 • WW 
IM1 = IMAX - 1 
I M2 - IMAX - 2 
KM1 = KMAX - 1 
KM2 *= KMAX - 2 
IF(ITN.EO.I) THEN 

C ! ! ! ! ! LOCAL TIME STEP OPTION FOR STEADY STATE CALCULATIONS 
IF(REDFRE.LT. 0.001) THEN 
DO 777 K - 2 . KMAX - 1 

DO 777 I = 2 . IMAX - 1 

DELTAT(I.K) = 0.5 / (1. + SORT (ABS(YACOB( I ,K) ) ) ) 

777 CONTINUE 
ELSE 

DO 778 K *= 2 , KMAX - 1 

DO 778 I = 2 . IMAX - 1 

778 DELTAT ( I ,K) = 1.0 
END IF 

END IF 
C 

C*** THE DISSIPATION TERMS ARE CONSTRUCTED AND STORED IN THE ARRAYS DQ1 . 

C»»» D02.DQ3 AND D04 . 

C 

CALL DISSIP 
C 

C*.* THE RIGHT HAND SIDE AT KNOWN TIME LEVEL IS NOW COMPUTED AND ADDED 
CALL RESI (OMEGA) 

C 

C*** IF VISCOUS FLOW IS COMPUTED THE VISCOUS TERMS ARE ADDED TO DQ1 ETC. HERE. 
C 

IFfREYREF.GT .0. ) CALL STRESS ( I TN, DALFA) 

C**. I— SWEEP . 

C 

DTH = DT • 0.5 
DTW = DT .WWI 
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DO 3 K - 2 , KM1 

CALL AMAT1 (K. [MAX-1 , XIX. XIZ.XIT) 

DO 4 II -1.4 

DO 4 12 -1.4 

DO 5 I - 2 . I MAX - 1 

EE( II .12. 1 — 1 . K ) - A( 1 1 , 12 , 1 + 1 ) • DTH • DELTAT(I.K) 
DD (11.12,1 — 1.K) - - A( 1 1 , 12 , 1 — 1 ) • DTH • DELTAT(I.K) 
5 CONTINUE 
4 CONTINUE 



C 

C*** IMPLICIT DAMPING ADDED HERE. 

C 



DO 6 II 
DO 6 I 
DTI 

DD( 1 1 .11. 1 — 1 . K) 
EE ( 1 1 .11,1-1 ,K) 
*44( I 1 . II . 1-1 ,K) 
6 CONTINUE 
DO 990 I 
GG( 1 . 1-1 .K) 
GG(2. 1-1 ,K) 

GG( 3 , 1-1 ,K) 

GG(4 . 1-1 ,K) 

990 CONTINUE 
3 CONTINUE 



1 . 4 

2 . I MAX — 1 

DTH / YACOB(I.K) • DELTAT ( I . K) 

DD( 1 1 . 1 1 , 1-1 .K) - DTI * YACOB(I-I.K) 
EE(I1,I1 .I-1.K) - DTI • YACOB(I+1.K) 
1 . + 2. • DTW . DELTAT( I ,K) 

2 , I MAX - 1 

DOI ( I ,K) • DELTAT(I ,K) 

D02 ( I ,K) • DELTAT ( I . K ) 

D03(I .K) • DELTAT ( I ,K) 

D04( I ,K) . DELTAT ( I ,K) 



C 

C 

C 
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PERFORM BLOCK TRIDIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 



CALL MATRX 1 ( I MAX . KMAX ) 



DO 991 K 
DO 991 I 
DQ1 ( I ,K) 
D02 ( I , K) 
DQ3( I ,K) 
D04( I , K) 
CONTINUE 



2 . KMAX - 
2 . I Ml 
GG( 1 . 1-1 .K) 
GG( 2 , 1-1 , K ) 
GG(3, 1-1 , K ) 
GG( 4 , 1-1 ,K) 
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C 

C K-SWEEP BEGINS HERE. 

C 

DO 13 I = 2 . I Ml 

CALL AMAT2 ( I . KMAX— 1 . ZETAX , ZETAZ . ZETAT) 

DO 15 II -1.4 

DO 15 12 —1.4 

DO 15 K = 2 . KMAX - 1 

EE( 1 1 . 1 2 . 1 . K— 1 ) - A(1 1 . I2.K+1 )*DTH • DELTAT ( I , K ) 
DD( 1 1 . 1 2 . I . K-1 ) - —A (11,12, K— 1 )*DTH • DELTAT ( I . K ) 
15 CONTINUE 



C 

C... SECOND ORDER DAMPING ADDED HERE. 

C 

DO 16 II 
DO 16 K 
DTI 

DD ( 1 1 .11.1 .K-1 ) 

EE ( 1 1 , 1 1 . 1 , K-1 ) 

16 MM( 11.11,1 , K— 1 ) 

DO 17 K 
GG ( 1 . I ,K— 1) 

GG(2. I , K— 1 ) 

GG(3 , I ,K-1 ) 

GG(4 , 1 , K— 1 ) 

17 CONTINUE 
13 CONTINUE 

C 

C PERFORM BLOCK TRI DIAGONAL MATRIX INVERSION FOR THE ENTIRE PLANE 

C 



-1.4 

- 2 . KMAX - 1 

= DTW / YACOB(I.K) • DELTAT(I.K) 

« DD( 1 1.1 1.1. K-1) - DTI • YACOB ( I , K— 1 ) 
« EE( 1 1.1 1.1. K-1) - DTI • YACOB ( I , K+1 ) 

- 1 . + 2. • DTYf • DELTAT( I ,K) 

= 2 . KMAX - 1 

«= DQI(I.K) 

-= D02 ( I . K) 

*= DQ3 ( I . K ) 

= D04( I ,K) 
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CALL MATRX2( IMAX.KMAX) 



DO 18 I 


= 2 , IMAX - 1 




DO 18 K 


= 2 , KM1 




DO 1 ( I . K ) 


= GG(1 . I . K— 1 ) 




DQ2( I ,K) 


= CG(2, I . K— 1 ) 




DQ3( I ,K) 


= GG ( 3 , I . K— 1 ) 




DQ4( I , K) 


- GG(4 , I , K-1 ) 




18 CONTINUE 






•• UPDATE FLOW 


VARIABLES AT INTERIOR POINTS. 


967 CONTINUE 






RMAX 


- 0 . 




RUMAX 


= 0 . 




RVMAX 


* 0 . 




EMAX 


= 0 . 




DO 995 K 


= 2 . KM1 




DO 19 I 


= 2 . IM1 




01(1. K) 


= 01 ( I ,K) + D01 ( I , K ) 


. YACOB( I .K) 


02(1. K) 


= 02 ( I .K) + D02 ( I , K ) 


• YACOB( I .K) 


03 ( I ,K) 


= Q3( I .K) + DQ3( I .K) 


• YACOB( I .K) 


04(1 ,K) 


= 04(1 ,K) + D04 ( I , K ) 


• Y ACOB ( I , K ) 



19 CONTINUE 



Cl!!!! DETERMINE WHERE IN FLOW FIELD DENSITY IS CHANGING MOST RAPIDLY 
DO 995 I = 2 . IMAX - 1 

IF (RMAX . LT . ABS(DQ1 ( I , K) • YACOB( I , K) ) ) THEN 
IR - I 

KR - K 

END IF 



RMAX = AMAX1 (RMAX ,ABS(D01 ( I .K) • YACOB(I.K))) 

RUMAX = AMAX 1 ( RUMAX . ABS ( D02 ( I . K ) • YACOB(I.K))) 

RVMAX = AMAX 1 ( RVMAX , ABS( DQ3( I . K) • YACOB(I.K))) 

EMAX = AMAX 1 ( EMAX , ABS ( D04 ( I . K ) • YACOB(I.K))) 

995 CONTINUE 

C I F ( ( I TN— 1 )/1 00* 1 00 . EQ . ( ITN— 1 )) WRITE (6.3002) 

C I F( ITN .EQ. 0) WRITE (6,3002) 

C ! ! ! ! ! SELECT INTERVAL AT WHICH OUTPUT OF RESIDUALS IS DESIRED 
C IF( ( ITN-1 )/10*10. EQ. ( ITN-1 ) ) WRITE (6,3001) RMAX , RUMAX , RVMAX , 

C 1EMAX.IR.KR 

RETURN 

3002 FORMAT (// . 4X . ' DRMAX ‘ . 1 1 X . * DUMAX ‘.11X,’ DVMAX ’.11X,’ DEMAX ’ . 1 0X . 

1 ’ IR’ . 3X . ’KR’ ) 

3001 FORMAT(4(E14.8.2X) ,215) 

END 



C 

C 

C 



C 

C 

C 

C 



C 



SUBROUTINE MATRX1 ( IMAX.KMAX) 

CCA440N/TR ID/DD( 4,4,161 . 41 ) ,MM(4 , 4 . 1 61 . 41 ) . EE(4. 4 . 1 61 .41 ) , 

1GG(4, 161 .41 ) 

COMMON /SCR AT /A( 4 , 4 , 161 ) ,HH(4 .4 , 1 61 . 41 ) ,C(4.5 . 161 ) 

R EAL MM 

REAL L1 1 . L21 . L31 . L41 . L22 . L32 . L42 , L33 . L43. L44 
2.L1I.L2I.L3I.L4I 

THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
AN ENTIRE PLANE DURING THE XI- SWEEP 



DO 1 11=1, 4 
DO 1 K = 2 , KMAX - 1 
A I = 1 . / MM( 1 , 1 , 1 . K) 
GG(I1 .1 ,K) = GG ( 1 1 , 1 , K ) 
HH( II ,1 . 1 ,K) 

HH( II ,2. 1 .K) 

HH( II ,3. 1 ,K) 

HH( II .4.1 . K ) 

1 CONTINUE 



A I 

EE( 1 1 .1 . 1 ,K) • Al 
EE (11.2.1 ,K) * AI 
EE( II .3. 1 ,K) • AI 
EE( II .4.1 , K ) • AI 
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DO 1000 I = 

DO 5 11- 

DO 2 K - 

C( 1 1 .1 .K) 

1 

2 

3 

2 CONTINUE 
DO 5 12 « 

DO 5 K - 

A ( 1 1 , 12, K)* 

1 

2 

3 



2 , I MAX - 2 
1 . 4 

o KUAX — 1 

- GG( II . I .K) - DD( 1 1 . 1 . I ,K) 

- DD( 1 1 .2. I , K) 

- DD( II ,3.1 , K ) 

- DD(I1 .4,1 ,K) 



1 , 4 

2 . KMAX - 1 

|yM(I1.I2.I.K) - DD(Il.l.l.K) 

- DD (II .2.1 ,K) 

- DD(I1 .3.1 ,K) 

- DD( 1 1 .4.1 .K) 



C( 1 1 . 12+1 .K)- EE( 1 1 . 12. 1 .K) 
5 CONTINUE 

DO 3 K - 2 . KMAX - 1 



• GG( 1 , 1-1 .K) 

• GG( 2 . 1-1 , K ) 

• GG( 3 . 1-1 .K) 

• GG(4. 1-1 , K ) 



• HH(1 , 12, 1-1 .K) 

• HH(2. 12. 1-1 , K ) 

• HH(3, 12. 1-1 ,K) 

• HH(4, 12. 1-1 .K) 



1 

1 



1 

1 



1 

1 



1 

1 



1 

1 



LI 1 = A( 1 , 1 , K ) 

LI I = 1 . / LI 1 
U 1 2 = A(1 ,2.K) • LI 1 
U13 «= A(1 ,3.K) • LI 1 
U1 4 = A(1 ,4.K) • LI 1 
L21 = A(2, 1 .K) 

L31 = A( 3 , 1 , K ) 

L41 = A( 4 , 1 .K) 

L22 - A( 2 . 2 , K ) - L21 • U12 
L2I = 1 . / L22 

U23 = ( A( 2 , 3 , K ) - L21 • U13) • L2I 

U24 = ( A( 2 . 4 . K) - L21 • U14) • L2I 

L32 = A(3. 2 , K) - L31 • U12 

L42 - A( 4 , 2 , K) - L41 • U12 

L33 - A(3 , 3 , K) - L31 • U13 - L32 • U23 

L3I - 1 . / L33 

U34 - ( A( 3 , 4 , K ) - L31 • U14 - L32 • U24) • L3I 

L43 - A( 4 , 3 , K ) - L41 • U13 - L42 • U23 

L44 = A(4 , 4 , K ) - L41 • U14 - L42 • U24 - L43 • U34 



L4 I — 1 . 


/ 


L44 




C(1 .1 .K) 




C( 1 . 1 , K ) 


• LI I 


C(2, 1 .K) 


= 


(C(2. 1 ,K) 


- L21 


C( 3, 1 , K) 


- 


(C(3. 1 , K ) 


- L31 








- L32 


C ( 4 1 .K) 




(C(4 , 1 . K) 


- L41 


c :.k) 


= 


C(1 ,2,K) 


• LI I 


C. .K) 


3= 


(C(2,2.K) 


- L2 1 


C. .K) 


*= 


(0(3. 2. K) 


— L31 








- L32 


C( 4 , 2 . K) 


* 


(C ( 4 , 2 . K ) 


- L41 


C(1 .3.K) 




C(1 . 3 . K) 


• LI 1 


C ( 2 , 3 . K) 


= 


(C(2,3,K) 


- L2 1 


C(3,3,K) 


3 


(C(3,3,K) 


— L31 








- L32 


C(4,3,K) 


* 


(C(4 ,3,K) 


- L41 


C(1 . 4 , K ) 


. 


C(1 . 4 . K ) 


• LI I 


C(2.4,K) 


« 


(C ( 2 . 4 . K ) 


- L2 1 


C(3.4,K) 




(C (3 , 4 , K ) 


— L31 








- L32 


C(4 , 4 , K) 


« 


(C( 4 . 4 , K ) 


- L41 



C(1 .1 .K)) • L2I 
C(1 .1 , K ) 

C(2 , 1 ,K) ) • L3I 

C(1 .1 ,K) - L42 • C(2 . 1 .K) 

- L43 • C(3 . 1 , K ) ) • L4I 

C(1 ,2,K)) • L2I 
C(1 .2.K) 

C(2 , 2 ,K) ) • L31 

C(1 , 2 ,K) - L42 • C( 2 . 2 , K) 

- L43 • C(3,2,K) ) • L4I 

C ( 1 ,3,K) ) • L2I 
C(1 .3.K) 

C(2 , 3. K) ) • L3I 

C(1 ,3,K) - L42 • C(2 , 3 . K) 

- L43 • C(3,3,K) ) • L4I 

C(1 ,4.K)) • L2I 
C(1 . 4 , K ) 

C( 2 . 4 , K ) ) • L3I 

C(1 , 4 , K ) - L42 • C(2 . 4 , K) 

- L43 • C(3 , 4 , K) ) * L4I 



C(1 . 5 . K) - 


C(1 ,5,K) 


• LI I 


C ( 2 , 5 , K ) = 


(C( 2 . 5 . K ) 


- L21 


C ( 3 . 5 , K ) = 


(C ( 3 . 5 , K ) 


- L31 






- L32 


C ( 4 . 5 . K) = 


(C( 4 , 5 , K ) 


- L41 



C(1 .5.K)) • L2I 
C(1 . 5 . K ) 

C (2 . 5 . K) ) • L3I 
C(1 . 5 . K) - L42 • C(2,5,K) 
- L43 • C(3 . 5 , K) ) • L4I 
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C f 3 . 1 ,K) 


- 


C(3 , 1 .K) 


- 


U34 


• 


C ( 4 , 1 , K ) 


C(2. 1 .K) 


* 


C ( 2 . 1 , K) 


- 


U24 


• 


C ( 4 , 1 , K) 


1 






- 


U23 


• 


C(3, 1 , K ) 


C(1 , 1 .K) 


« 


C(1 .I.K) 


- 


U14 


* 


C ( 4 , 1 , K) 


1 






- 


U13 


* 


C ( 3 . 1 , K ) 


C(3 . 2 . K 1 




C(3 , 2 . K) 


- 


U3* 


• 


C( 4 , 2 , K) 


C(2.2,K) 


« 


C(2.2,K) 


- 


U24 


• 


C( 4 , 2 , K) 


1 






- 


U23 


• 


C( 3 , 2 . K) 


C(1 , 2 . K) 


■ 


C(1 .2.K) 


- 


U14 


• 


C( 4 , 2 . K ) 


1 






- 


U13 


• 


C( 3 , 2 , K) 


C(3 , 3 . K) 


- 


C(3 , 3 . K) 


- 


U34 


• 


C( 4 , 3 , K) 


C( 2 . 3 , K) 




C(2.3.K) 


- 


U24 


• 


C(4 , 3 , K) 


1 






- 


U23 


• 


C(3.3.K) 


C(1 .3.K) 




C( 1 ,3,K) 




U14 


« 


C ( 4 . 3 . K) 


1 






- 


U13 


• 


C( 3 . 3 . K) 


C ( 3 , 4 , K ) 


- 


C(3.4,K) 


- 


U34 


• 


C ( 4 , 4 , K) 


C ( 2 . 4 , K ) 


■ 


C(2 , 4 , K) 


- 


U24 


• 


C ( 4 , 4 , K) 


1 






- 


U23 


• 


C( 3 , 4 . K) 


C(1 .4.K) 




C(1 . 4 , K) 


- 


U14 


• 


C ( 4 . 4 , K ) 


1 






- 


U13 


* 


C ( 3 , 4 . K ) 


C(3 . 5 . K) 


- 


C(3 . 5 . K) 


- 


U34 


m 


C ( 4 , 5 , K ) 


C(2 , 5 , K) 


« 


C(2 , 5 , K) 


- 


U24 


* 


C( 4 . 5 . K) 


1 






- 


U23 


* 


C( 3 , 5 , K) 


C(1 , 3 ,K) 


« 


C ( 1 . 5 , K) 


- 


U14 


m 


C ( 4 , 5 , K ) 


1 






- 


U13 


m 


C( 3 . 5 . K) 


3 CONTINUE 














DO 6 11 




= 1 . 


4 








DO 9 K 




- 2 . 


KMAX - 


■ 1 




9 GG( 11 . I ,K) 


- C( 1 1 


.I.K) 






DO 6 12 




- 1 . 


4 








DO 6 K 




= 2 . 


KMAX - 


• 1 




HH( 11 . 12. 


I.K) - C( I 1 


. 12+1 . 


K) 





6 CONTINUE 
1000 CONTINUE 



C 

C 

C 



BACKWARD SUBSTITUTION 



DO 7 I 
DO 7 II 
DO 7 K 
GG( 1 1 , I . K) 

1 

2 

3 

7 CONTINUE 
RETURN 
END 



I MAX -3,1 . - 1 

1 . 4 

2 . KMAX - 1 

GG(II.I.K) - HH( 1 1 , 1 , I ,K) 

- HH( 1 1 ,2, I ,K) 

- HH( 1 1 ,3, I , K ) 

- HH( 1 1 ,4, I , K) 



U12 • C(2. 1 , K ) 



U12 • C(2 . 2 . K) 



U12 • C(2 , 3 , K) 



U12 • C(2, 4 ,K) 



U12 • C(2 . 5 , K) 



• GG( 1 . 1+1 , K) 

• GG(2. 1+1 , K) 

• GG(3. 1+1 , K) 

• GG(4. 1+1 , K) 



C 

C 

C 



C 

C 

C 

C 



SUBROUTINE MATRX2( IMAX .KMAX) 

C0fc*ADN/TR I D/DD (4,4, 161 . 41 ) ,MM(4 . 4, 1 61 .41 ) ,EE(4,4, 161 .41) , 

1GG(4. 161 .41 ) 

C0k*40N/SCRAT/A(4,4.161) ,HH(4.4,161 . 41 ) ,C(4 ,5. 161 ) 

REAL MM 

REAL L1 1 . L21 . L31 . L41 . L22 . L32 . L42 . L33 . L43 . L44 
2.L1I.L2I .L3I.L4I 

THIS SUBROUTINE PERFORMS THE BLOCK TRIDIAGONAL MATRIX IVERSION FOR 
AN ENTIRE J-CONSTANT PLANE DURING THE 2 ETA- SWEEP 



DO 1 II 
DO 1 I 
AI 

GG( 11.1,1) 
HH( 11.1,1.1) 



1 . 4 

2 . IMAX - 1 

1 . / 1 * 1 ( 1 . 1 , 1 , 1 ) 
GG (II .1.1) • AI 
EE(I1 .1 .1.1) • AI 
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HH( II .2. I . 1 ) 
HH( I 1 .3, I . 1 ) 
HH( I 1 .4.1,1) 
1 CONTINUE 



EEC 1 1 .2.1.1) • AI 
EE ( 1 1 .3,1,1) • AI 
EE( 1 1 .4. 1 . 1) • AI 



00 1000 K - 2 , KMAX - 2 

00 5 11-1,4 

00 2 I - 2 , I MAX - 1 

C(I1.1.I) - GG(II.I.K) - 00(1 1 , 1 . I ,K) • GG(1,I,K-1) 

1 - 00(11 , 2 , 1 ,K) • GG ( 2 , 1 , K— 1 ) 

2 - 00(11 , 3 , I ,K) • GG(3.I.K-1) 

3 - 00(1 1 ,4, I ,K) • GG ( 4 . I ,K-1 ) 

2 CONTINUE 

DO 5 12 -= 1 .4 

DO 5 I - 2 , IMAX - 1 

AC 1 1.12. 1)- WW(H .I2.I.K) - 00(1 1,1,1, K) • HH(1 ,12, I ,K-1 ) 

1 - 00(11, 2, 1. K) • HH(2, 12 , I ,K-1 ) 

2 - 00(11, 3. I, K) • HH(3, 12, I ,K-1 ) 

3 - 00(11, 4. 1. K) • HH(4. 12, I ,K-1 ) 

C( 1 1 . 12+1 , I)= EE( II . 12, I , K) 

5 CONTINUE 

DO 3 I - 2 . IMAX - 1 
L1 1 = A(1 .1 ,1) 

LI I = 1 . / L1 1 
U12 - A(1 ,2, I) 



LI I 
LI I 
LI I 



U 13 = A( 1 , 3 , 1 ) 

U 1 4 = A( 1 . 4 , I) 

L 21 -= A( 2. 1 , I) 

L 31 ■= A( 3 . 1 . 1 ) 

L 41 = A( 4 . 1 . I) 

L 22 - A( 2 . 2 . I) - L 21 • U 12 
L 2 I - 1 . / L 22 

U 23 «= ( A( 2 , 3 . 1 ) - L 21 • U 13 ) • L 2 I 

U 24 - ( A( 2 . 4 , I ) - L 21 • U 14 ) . L 2 I 

L 32 - A( 3 . 2 , I ) - L 31 • U 12 

L 42 ■= A( 4 , 2 . I) - L 41 • U 12 

L 33 - A( 3 , 3 , I) - L 31 • U 13 - L 32 • 

L 3 I - 1 . / L 33 

U 34 - ( A( 3 , 4 , I ) - L 31 • U 14 - L 32 
L 43 = A( 4 , 3 , I ) - L 41 • U 13 - L 42 
L 44 - A( 4 , 4 , I ) - L 41 • U 14 - L 42 



U23 



L3I 



• U24) 

• U23 

• U24 - L43 • U34 



L4I - 1 . 


/ 


L44 






C( 1 . 1 . I ) 




0(1 .1 . I) 


• LI I 




C (2 . 1 . I) 


= 


(0(2. 1.1) 


- L21 


• 


C(3 , 1 , I) 


= 


(0(3,1 . I) 


- L31 


• 


! 






- L32 


• 


C ( 4 , 1 , I ) 


* 


(0(4.1 . I) 


- L41 


* 


c (1 . 2 , 1 ) 


SE 


0(1 .2.1) 


• LI I 




C(2,2, I) 


= 


(0(2.2. I) 


- L2 1 


• 


C(3 . 2 . I ) 


= 


(0(3.2. I) 


— L31 


* 








- L32 


• 


C(4.2. I) 


= 


(0(4.2. I) 


- L41 


* 


C(1 .3.1) 


= 


0(1 ,3. I) 


• LI I 




C(2,3, I) 


= 


(0(2.3. I) 


- L21 


• 


C(3,3. I) 




(0(3.3, I ) 


- L31 


• 








- L32 


• 


C(4 , 3 , I ) 




(0(4.3. I) 


- L41 


• 


C(1.4, I) 




0(1 .4.1) 


• LI I 




0(2.4. I) 


= 


(0(2.4. I) 


- L21 


* 


C(3 , 4 , I ) 




(0(3.4. I) 


- L31 


* 








- L32 


• 


0(4,4, I ) 


- 


(0(4.4. I) 


- L41 


* 


0(1 .5.1) 


« 


0(1.5. I) 


• LI I 





C(1.1.I)) 
C(1 .1 .1) 
C( 2 , 1 , I ) ) 



L2I 



L3I 



- L43 • C(3 , 1 . 1)) • L4I 



C(1 .2.1) - L42 • C(2.2, I) 
- L43 * C(3 , 2 , I )) • L4I 



C(1 .3.1)) • 
C( 1,3,1) 

C(2 ,3. I)) • 



L2I 



L3I 



- L43 . C(3, 3,1)) • L4I 

C( 1 .4.1)) . L2I 
C ( 1 , 4 , I ) 

C ( 2 , 4 , I )) • L3I 

C(1 .4.1) - L42 • C(2 ,4,1) 

- L43 • C(3 , 4 , I )) . L4I 
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C(2 , 5 . I ) 


- (C( 2 , 5 . I 


) - L21 


C(3 . 5 , I ) 


- (C( 3 , 5 , I , 


) - L31 


1 




- 132 


C(4 , 5 . 1 ) 


- (C(4 , 5 , I ) - L41 


C(3 , 1 .1) 


- C( 3 , 1 .1) 


- U34 


C(2 , 1 .1) 


- C( 2 , 1 , I ) 


- U24 


1 




- U23 


C( 1 , 1 , 1 ) 


- C(1 .1.1) 


- U14 


1 




- U13 


C(3.2, I) 


- 0(3,2. I) 


- U34 


C(2.2. I) 


- C(2.2. I) 


- U24 


1 




- U23 


C( 1 , 2 , I ) 


= C(1 .2.1) 


- U 1 4 


1 




- U13 


C(3 , 3 , I ) 


- C ( 3 , 3 , I) 


- U34 


C(2.3. I) 


= C(2.3. I) 


- U24 


1 




- U23 


C(1 .3, I) 


- C(1 ,3. I) 


- U 1 4 


1 




- U13 


C(3 , 4 , I) 


- C(3.4. I) 


- U34 


C ( 2 . 4 , I ) 


- C (2 , 4 , I ) 


- U24 


1 




- U23 


C( 1 . 4 , I ) 


- C(1.4. 1) 


- U 1 4 


1 




- U13 


C(3 , 5 , I ) 


- 0(3.5, I ) 


- U34 


0(2,5, I) 


= 0(2,5, I) 


- U24 


1 




- U23 


C(1 .5, I) 


= C(1 .5. I) 


- U 1 4 


1 




- U 1 3 


3 CONTINUE 






DO 6 11 


- 1 . 


4 


DO 9 I 


= 2 . 


I MAX - 


9 GG( 1 1 . I , K) «= C ( 1 1 


.1.1) 


DO 6 12 


= 1 , 


4 


DO 6 I 


- 2 , 


IMAX - 


HH(I1 , 12. 


I . K) - C ( I 1 


.12+1.1 



6 CONTINUE 
1000 CONTINUE 



• C ( 1 .5.1)) • L2I 

• C(1 .5. I) 

• C(2.5. I)) • 13 1 

• C(1 .5. I) - L42 • C(2 .5 , I ) 

- 143 • C(3,5, I)) • L4I 
C(4, 1 , I ) 

0(4, 1 ,1) 

C(3.1 ,1) 

C(4 , 1 , I ) 

C ( 3 . 1 ,1) - U12 • C( 2 , 1 , I ) 
0(4,2, I ) 

C(4,2. I) 

C(3.2 , I ) 

C ( 4 , 2 . 1 ) 

C ( 3 . 2 . I ) - U12 • C(2 ,2 , I ) 

C( 4 , 3 , I ) 

C( 4 , 3 , I ) 

C(3.3. I) 

C(4.3. I) 

C ( 3 , 3 , I ) - U12 • C ( 2 , 3 , I ) 

C( 4 , 4 , I ) 

C ( 4 , 4 , I ) 

C( 3 , 4 , I ) 

C ( 4 , 4 . I ) 

C (3 , 4 , I ) - U12 • C(2 , 4 , I) 
C(4.5. I) 

C(4,5, I ) 

C( 3 , 5 , I ) 

C( 4 , 5 , I ) 

C(3 , 5 , I ) - U12 * C(2,5, I ) 



BACKWARD SUBSTITUTION 

DO 7 K = KMAX - 3, 1 . - 1 

DO 7 II -1,4 

DO 7 I - 2 , I MAX - 1 

GG(II.I.K) - GG(II.I.K) - m(I 1 , 1 ,I,K) • GG(1,1,K41) 

1 - HH( 1 1 ,2, I ,K) • GG ( 2 , I , K+1 ) 

2 - HH( 1 1 .3. I ,K) • GG(3 , 1 . K+1 ) 

3 - 1*1(11 ,4,1 ,K) • GG( 4,1, K+1) 

7 CONTINUE 

RETURN 

END 



SUBROUTINE METRIC 
C0A4ADN/F I X/OMEGA 

COMMON/DGRID/DT, I MAX , KMAX , ITEL . I TEU 
COMMON/GR I D 1 /X ( 1 6 1 , 41 ) , Z(161 ,41 ) 

C0A4ADN/GR I D/YACOB( 1 61 .41) 

COMMON/MTR I X/X I X ( 1 6 1 ,41 ) ,XI2( 161 .41 ) ,ZETAX(161 .41 ), 

+ZETAZ( 161 ,41 ) . 

1 XIT( 1 61 ,41) , ZETAT (161 ,41) 

•• SUBROUTINE METRIC COMPUTES THE METRICS IN BOTH DIRECTIONS AND 
THE UNSTEADY COEFFICIENTS ETAT . ETC. 
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c 

DO 1000 K - 1 , KMAX 
DO 1000 I - 1 . I MAX 
XTAU - OMEGA . Z(I ,K) 

YTAU = OMEGA . (-X ( I . K ) ) 

C • • • PRESENT SET UP IS TOR PLOW PAST AN AIRPOIL. 

C 

C ! ! ! ! (CENTRAL DIFFERENCES AT INTERIOR POINTS. TWO-POINT ONE-S.DED 
C! ! ! ! IDIFFERENCES DOWNSTREAM, THREE-POINT AT OTHER OUTER BOUNDARIES 
IF( I . EO. 1 .OR. I . EQ. IMAX) GO TO 10 
XXI * .5 • (X(I+1 , K)— X( 1-1 , K ) ) 

ZXI = .5 • (Z(I+1 . K)-Z( 1-1 ,K)) 

GO TO 15 

10 I F ( I .EO. IMAX) GO TO 16 

XXI « 1.0 • (X ( 2 , K ) - X(1 , K ) ) 

ZXI = 1.0* ( Z ( 2 , K ) - Z( 1 , K) ) 

GO TO 15 

16 XXI = 1.0 • (X(IMAX.K) - X( IMAX— 1 .K) ) 

ZXI = 1.0 • (Z(IMAX.K) - Z ( IMAX— 1 , K) ) 

15 CONTINUE 

IF(K. EO. 1 .OR.K. EQ.KMAX) GO TO 1 7 
XZET = .5 • ( X( I , K+1 ) — X ( I . K-1 ) ) 

ZZET = .5 • ( Z ( I , K+ 1 )— Z( I ,K— 1 )) 

GO TO 20 

17 IF(K. EQ.KMAX) GO TO 18 

XZET •= 2. • X( I . 2 ) — 1 .5 • X ( I . 1 ) - .5 • X(I,3) 

ZZET = 2. • Z( I .2) - 1.5 • Z ( I . 1 ) - .5 • Z(I,3) 

GO TO 20 

18 XZET - 1.5 • X( I , KMAX)— 2 . • X( I .KMAX-1 )+.5*X( I .KMAX-2) 

ZZET * 1.5 • Z ( I , KMAX)— 2 . • Z ( I , KMAX-1 )+. 5*Z ( I . KMAX-2 ) 

20 CONTINUE 

C! ! ! ! ! COMPUTE JACOBIAN 

YACOBI = XXI • ZZET - XZET • ZXI 
YACOB( I .K) = 1 . / YACOBI 
XIX(I.K) = ZZET • YACOB(I.K) 

XIZ(I.K) = -XZET • YACOB(I.K) 

XTAU - OMEGA . Z(I ,K) 

YTAU = - OMEGA • X ( I ,K) 

XIT(I.K) - - XIX(I.K) • XTAU - XIZ(I.K) • YTAU 
ZETAX(I.K) *= -ZXI » YACOB(I.K) 

ZETAZ(I.K) - XXI » YACOB(I.K) 

ZETAT(I.K) -= - ZETAX(I.K) • XTAU - ZETAZ(I.K) • YTAU 
1000 CONTINUE 
RETURN 
END 



SUBROUTINE DISSIP 

COMMON/ FLOW/Q1 ( 1 61 . 41 ) ,Q2( 1 61 . 41 ) ,Q3( 1 61 . 41 ) . Q4( 1 6 1 . 41 ) 
C0M40N/PERTR/DQ1 (161 , 41 ) , DQ2 ( 1 61 . 41 ) .003(161.41 ) ,DQ4( 161 , 41 ) 
COMMON/DGR I D/DT . IMAX . KMAX . I TEL . ITEU 
C0U4ON/GR I D/YACOB (161.41) 

C0M4ON/DAMP/WW . WW I 

DIMENSION P( 1 61 ) , EPS( 161 ) . DIS1 ( 1 61 . 4) . DIS2 ( 1 61 , 4) 

THIS SUBROUTINE ADDS THE FOURTH ORDER DISSIPATION TERMS TO THE 
RIGHT HAND SIDE 

I Ml = IMAX - 1 
KM1 = KMAX - 1 
I M2 = IMAX - 2 
KM2 = KMAX - 2 

C 

DO 10 K=2 . KM1 

C COMPUTE SWTICHING FUNCTION BASED ON SECOND DERIVATIVE OF PRESSURE 
DO 1 I = 1 . IMAX 
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1 P(I) - .4 • (04(1. K)-(Q2(I ,K)«. 2+03(1, K) • «2)/(2 . *01 ( I . K ) ) ) 
DO 2 I -1 . I MAX 

IP2 - I + 2 

I F( I .EQ. IM1) IP2 - IMAX 
IM2 - I - 2 
I F ( I . EQ . 2 ) I M2 « 1 
I P 1 “1 + 1 

I F ( I .EQ. IMAX) IP1 - IMAX 
IM - I - 1 
I F( I .EQ. 1 ) IM - 1 
IF(I.EQ.I) I M2 - 1 
I F( I .EQ. IMAX) IP2 - IMAX 

EPS(I) - ABS(P(IP1)-2.«P(I)+P(IM))/ABS(P(IP1)+2.«P(I)+P(IM)) 
C VOL - 2. /(YACOB( I , K)+YACOB( I Pi ,K) ) 

VOL = 1 . 

DISI(I.I) - (01 ( I P 1 , K )— Q1 (I . K) ) »VOL 

DIS1 ( I .2) - (Q2(IP1 ,K)-Q2(I ,K)).VOL 

D1S1 ( I . 3) - (Q3(IPl .K)-03(I ,K)).VOL 
HP 1 = 04 ( I PI , K)+P( I Pi ) 

HP - 04 ( I , K)+P( I ) 

HM1 = Q4( IM.K) + P( IM) 

HP2 = Q4( IP2.K) + P( IP2) 

HPM = 04 ( IM,K)+P( IM) 

DIS1 ( I , 4) = ( HP 1 —HP ) • VOL 

0IS2 ( I . 1 ) =— (01(1 P2 , K)-3 . • (Q1 ( I PI ,K)-01(I.K))-01( IM.K) ) .VOL 
DIS2( I .2) — (Q2( I P2 . K)-3 . • (Q2( I PI ,K)— Q2( I ,K) )— Q2( IM . K) ) »VOL 
DIS20.3) =- (Q3( IP2,K)-3. *(Q3( I PI ,K)— Q3( I , K))-Q3( IM,K))*VOL 
DIS2 ( I .4) — ( HP2-3 . • ( HP 1 —HP)— HPM) • VOL 
2 CONTINUE 

DO 15 I * 1 . IM1 

D2P = AMAX1 ( EPS( I ) . EPS( 1+1 )) 

C22 = 60. • D2P 



C11 = AMAX1 (0.0. ( 1 .-C22) ) 

C22 = C22 • WW/YACOB( I ,K) . DT 
C11 - C11 • WW/YACOB( I ,K) • DT 

C ! ! ! ! ! SW1 TCH ON SECOND-ORDER AND SWITCH OFF FOURTH-ORDER DISSIPATION 
C ! ! ! ! ! IN VICINITY OF SHOCKS 

DISI(I.I) = C11 • DIS2( I , 1 ) + C22 * DISI(I.I) 

DIS1 ( I .2) - C11 • DIS2( I ,2) + C22 • DIS1(I,2) 

DIS1 ( I .3) = C11 • D I S2 ( I ,3) + C22 • DIS1(I,3) 

DIS1 ( I .4) = C11 • DIS2( I ,4) + C22 • DIS1(I.4) 

15 CONTINUE 

DO 16 I = 2 . IM1 

DQI(I.K) = DISI(I.I) - DIS1 ( 1-1 , 1 ) 

DQ2 ( I . K ) = DIS1 (1,2) - DIS1 ( 1—1 ,2) 

DQ3 ( I . K ) - DIS1 ( I ,3) - DIS1 ( 1—1 ,3) 

DQ4(I.K) - DI SI ( I . 4) - DI SI ( 1-1 , 4) 

16 CONTINUE 
10 CONTINUE 

C K DIRECTION 

C! !!! .'FOURTH-ORDER DISSIPATION ONLY 
DO 30 I - 2 . I M 1 
WT= 0.5 • DT • WW / YAC0B(I,2) 

W3 « 0.5 * DT . WW / YACOB(I.KMI) 

DQ1 (1,2) -WT. (01(1,1) - 2. • 01(1.2) + 01(1,3)) 

1+D01 ( I .2) 

002(1.2) -WT* (02(1.1) - 2. . 02(1.2) + 02(1.3)) 

1 +DQ2 ( I .2) 

DQ3( 1.2) -WT. (03(1.1) - 2. • 03(1.2) + Q3(1.3)) 

1 +D03 ( I .2) 

004(1,2) -WT. (04(1.1) - 2. • 04(1.2) + 04(1.3)) 

1+DQ4 ( I .2) 

WT= W3 



DQI(I.KMI) -WT. (01(1. KM2) - 2. • 01(1. KM1) + 01(1, KMAX ) ) 
1+DQ1 (I , KM1 ) 

DQ2( I , KM1 ) =WT • (02( I , KM2) -2. * Q2(I.KM1) +02(1, KMAX)) 
1 +D02 ( I . KM1 ) 
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DQ3( I . KM1 ) -WT. (Q3( I . KM2 ) - 2. • Q3(I.KM1) + Q3(I.KMAX)) 
1+DQ3( I . KM1 ) 

D04(I.KM1) «WT» (04(1, KM2) - 2. • 04(1, KM1 ) + 04(1, KMAX)) 
1+DQ4(I , KM1 ) 

00 35 K - 3 , KW2 

WT- - 0T • WW / YAC06( I . K ) 

OQI(I.K) -WT. (01(1, K+2) - 4. . 01(1, K+1) + 6. • Q1( 

II. K) - 4. • 01(1. K-1) + 01(1, K-2) )+001 ( I , K ) 

002(1. K) -WT*(Q2( I ,K+2) - 4. . 02(1. K+1) + 6. • Q2( 

II. K) - 4. • 02(1, K-1) + 02(1, K-2 ) )+DQ2 ( I . K ) 

003(1- K) -WT *(03(I . K+2) - 4. • Q3(I,K+1) + 6. • Q3( 

II, K) - 4. . 03(1, K-1) + Q3( I , K-2 ) )+DQ3( I , K ) 

004(1, K) -WT.(04( I ,K+2) - 4. • 04(1, K+1) + 6. • Q4( 

II. K) - 4. . 04(1. K-1) + 04 ( I ,K-2))+DQ4(I ,K) 

35 CONTINUE 
30 CONTINUE 

RETURN 

END 



SUBROUTINE WALLBC 
COM*DN/SURF/PSUR( 1 61 ) 

C0M4ON/GRID1/X(161 ,41).Z(161.41) 

COM40N/PAR/GAU4A . REYREF . ALFA . ALFA1 . REOFRE . AMI NF . ALFA I 
C0M40N/DGRID/DT. IMAX.KMAX, ITEL, ITEU 
C0fc#4ON/GR I D/YAC06( 1 61 ,41) 

CCM40N/DAMP/WW , WW I 

C0»44ON/MTRIX/XIX( 161 , 41 ) . X I Z ( 1 61 , 41 ) ,ZETAX( 161 .41 ) , 

+ZETAZ( 1 61 .41 ), 

1 X IT ( 1 61 .41 ) , ZETAT ( 1 61 ,41) 

C0fc#4ON/FLOW/Q1 ( 1 61 , 41 ) , 02 ( 1 61 , 41 ) , Q3( 1 61 , 41 ) . Q4( 1 61 . 41 ) 
DIMENSION Cl (4) 

DIMENSION A(2,2) ,RHS(2) 

C ! ! ! ! ! APPLY BOUNDARY CONDITIONS ON THE CUT AND THE AIRFOIL SURFACE 
DO 9 I-ITEU. IMAX 
II = IMAX +1-1 

01(1.1) = .5 • (01(1.2)401(11.2)) 

02(1,1) -.5* (02(1 .2)402(11 .2)) 

03(1,1) - .5 • (03(1,2) + 03(11,2)) 

04(1,1) = .5 • (04(1,2)404(11,2)) 

01(11 , 1 )=Q1 (1.1) 

02(11,1 )«Q2( 1.1) 

03(11 , 1 )=03( 1,1) 

04(11 , 1 )=04( 1.1) 

9 CONTINUE 

DO 1 1= ITEL . ITEU 
K = 3 

C 1 ( 1 ) - XIT(I.K) 

Cl ( 2) = XIX(I.K) 

Cl (3) = XIZ(I.K) 

UC0N3 = (Q2( I ,K)*C1 (2)+03( I ,K)*C1 (3)) 

1/01 ( I , K ) 

K-2 

Cl ( 1 ) = XIT(I.K) 

Cl (2) = XIX( I ,K ) 

Cl ( 3) - XIZ(I.K) 

UC0N2 - (Q2(I.K)-C1 ( 2)+Q3( I ,K)*C1(3)) 

1/01(1 , K ) 

RHS( 1 ) = 2. * UC0N2 - UC0N3 - XIT(I.I) 

C FOR VISCOUS FLOWS SET UCON TO ZERO ALSO 

I F(REYREF . GT .0. ) RHS(1) - - XIT(I.I) 

A ( 1 , 1 ) - X I X( 1 . 1 ) 

A(1 .2) = X IZ( I .1) 

A(2.1) = ZETAX(I.I) 

A(2 . 2) = ZETAZ(I.I) 
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RHS ( 2 ) - - ZETAT (1.1) 

TEMPI - A ( 1 , 1 ) 

TEMP2 - A(1 .2) 

TEMP3 - A(2. 1) 

TEMP* - A(2 , 2) 

DEN - 1. /(TEMPI • TEMP 4 - TEMP2 • TEMP3) 

A(1 .1) - A(2 . 2) • DEN 

A(1 .2) - - TEMP2 • DEN 

A(2, 1 ) - - TEMP3 • DEN 

A(2 . 2) - TEMPI • DEN 

01(1.1) - 2. • 01(1.2) - 01(1.3) 

02(1.1) = 01 ( I . 1 ) • (A( 1 , 1 ).RHS( 1 )+A(1 ,2)*RHS(2)) 

03(1.1) -01(1.1) • (A(2 . 1 ) »RHS( 1 )+A( 2 . 2) *RHS(2) ) 

1 CONTINUE 

DO 10 I = I T EL . I TEU 
U2=02( 1.2) /O 1(1 .2) 

*V2=03(I.2)/01(I .2) 

P2=(GAK*1A-1 . ).(Q4( I ,2)-0.5»O1 ( 1.2) • (U2»U2+W2»W2) ) 

U3=02 ( I . 3 )/Q1 ( I .3) 

W3=Q3( I . 3)/Q1 (1.3) 

P3=(GAMVA-1 . ) • (Q4( I , 3)-0 . 5*01 ( I , 3) • (U3»U3+W3*W3) ) 

P1*=( 4 . »P2-P3)/3 . 

PSUR ( I )-(GAA*4A.P1-1 . )/( . 7*AMINF**2) 

Ul=02( I . 1 )/Q1 (1.1) 

V»1=Q3( I . 1 )/Q1 (1.1) 

10 04(1 . 1 )-P1/(GAMMA-1 . )+0 . 5*01 ( I . 1 ) • (U1 -U1+W1 *W1 ) 

RETURN 

END 

C 



c 

SUBROUTINE STRESS( ITN.DALFA) 

OXMON/FLOW/Q1 (161 .41). 02(161 .41) ,03 (’61 .41) .04(161 .41 ) 
CC4440N/DGR I D/DT . I MAX . KMAX . I TEL. I TEU 
CC4440N/GR I D1 /X (161 .41),Z(161.41) 

COMvtON/PAR/GAUM.REYREF.ALFA.ALFAl , REDFRE . AMINF . ALFA I 
C0*440N/PERTR/D01 (161 . 41 ) , D02 ( 1 61 , 41 ) .003(161 .41 ) ,DQ4( 1 61 .41 ) 
C0*440N/MU TUR/CMU (161,41) 

DIMENSION AA( 1 61 .41 ) 

1 .RH1 (161 ) .RH2( 161 ) ,RH3(161 ) .RH4( 161) 

C0*440N/ LOG I C/RSTRT . P I TCH . P LUNG E 
LOGICAL RSTRT. PI TCH. PLUNGE 
U(I.J) = 02(1. J) / 01(1, J) 

V(I.J) = 03 ( I . J ) / 01(1. J) 

C THIS SUBROUTINE ADDS VISCOUS TERMS TO THE RIGHT HAND SIDE 
GOGM = GAMMA / (GAH4A - 1 . ) 

IF(ITN.GT. 10. OR. (RSTRT)) CALL EDDY(DALFA) 

C COMPUTE U AND V 

KMAXM1 = KMAX - 1 
IMAXM1 = IMAX - 1 
PR = 1 . 

DO 10 K « 1 . KMAX 

DO 10 I « 1 . IMAX 

E - Q4( I , K ) / 01(1. K) - 0.5 • (U( I , K) • *2+V( I . K) • *2) 

10 AA(I ,K) - GOGM . E 
C 

C COMPUTE TXX.TXY AND VISCOUS DISSIPATION AT I - 1 / 2 
C 

DO 30 K - 2 . KMAXM1 

DO 20 I - 2 . IMAX 

UXI - U(I , K) - U ( I — 1 .K) 

VXI = V( I .K) - V( 1-1 , K) 

AX I = AA( I ,K) - AA( 1—1 .K) 

UZET= . 25» (U( I , K+1 )— U( I ,K-1)+U( 1-1 .K+1 )— U( 1—1 ,K-1)) 

VZET= . 25* (V( I , K+1 )— V( I ,K-1 )+V( 1-1 .K+1 )-V( 1-1 .K-1 ) ) 

AZET= ,25»(AA( I .K+1 )-AA( I ,K-1 )+AA(I-1 .K+1 )-AA( 1-1 .K-1 )) 

XXI - X(I .K) - X ( I— 1 .K) 
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ZXI = Z( I .K) - Z( 1-1 , K ) 

XZET- .25 • ( X( 1 , K+1 )-X( I , K— 1 )+X ( I— 1 , K+ 1 )— X ( 1-1 , K- 1 ) ) 
ZZET- .25 • ( Z( I , K+1 )— Z( I . K— 1 )+Z( I — 1 , K+1 )— Z( I — 1 . K — 1 ) ) 

YAC - XXI • ZZET - ZXI • XZET 
YAC - 1 . / YAC 
XIX - ZZET • YAC 
ZETAX- - ZXI • YAC 
XIZ = -XZET • YAC 
ZETAZ- XXI • YAC 

CNM = .5 • (CMU(I.K) + CMU(I-I.K)) 

UX = UXI • XIX + UZET • ZETAX 
VX * VXI • XIX + VZET . ZETAX 
AX = AX I • XIX + AZET • ZETAX 

UZ = UXI • XIZ + UZET . ZETAZ 

VZ = VXI • XIZ + VZET • ZETAZ 

AZ - AX I • XIZ + AZET • ZETAZ 

TXX - —(—4 . • UX + 2. • VZ) • CNM / 3. 

TXY = CNM . (UZ + VX) 

TYY - -CNM / 3. • (-4. • VZ + 2. • UX) 

R4 = ( ( U ( I . K ) +U (1 — 1 , K ) ) *TXX+( V( I . K)+V( I— 1 ,K)) »TXY) »0 . 5 
1 + CNM / PR/(GAK*M - 1.) • AX 

S4 = ((U(I . K )+U( 1—1 .K))*TXY+(V(I .K)+V(I-1 ,K) ).TYY)*0 . 5 
1 + CNM / PR / (GAMMA - 1 . ) • AZ 

DEBUG 

TURN OFF ENRGY DISSIPATION AND DIFFUSION 
R4 - 0. 

S4 = 0. 

RH1 ( I ) = 0. 

RH2 ( I ) = (XIX . TXX + XIZ • TXY) / YAC 
RH3 ( I ) - (XIX . TXY + XIZ • TYY) / YAC 
20 RH4 ( I ) = (XIX • R4 + XIZ • S4) / YAC 
DO 30 I = 2 , IMAXM1 

DQI(I.K) = DO 1(1 , K ) + RH1 ( 1+1 ) - RH1 ( I ) 

D02 (I , K ) = D02( I , K) + RH2(I+1) - RH2( I ) 

D03 ( I , K ) - DQ3( I , K ) + RH3(I+1) - RH3(I) 

D04(I.K) = D04 ( I , K ) + RH4( 1+1 ) - RH4(I) 

30 CONTINUE 
C IN THE Z DIRECTION 

DO 70 I = 2 . IMAXM1 

DO 60 K = 2 . KMAX 

UXI = .25 • (U( 1 + 1 ,K)-U( 1 — 1 ,K )+U( 1 + 1 , K— 1 )— U( I— 1 , K— 1 ) ) 

VXI = .25 • ( V( 1 + 1 , K )— V ( I — 1 , K )+V ( 1 + 1 , K— 1 )— V ( I — 1 , K— 1 ) ) 

AX I = .25 • ( AA( 1 + 1 , K )— AA( I — 1 . K )+AA( 1+1 , K— 1 )— AA( I— 1 . K— 1 ) ) 
XXI = .25 • (X(I + 1,K)-X(I — 1.K )+X ( 1 + 1 , K— 1 )— X ( 1 — 1 ,K— 1 ) ) 

ZXI = .25 • (Z( 1 + 1 . K ) — Z (I — 1 ,K)+Z(I + 1,K— 1 )— Z( 1-1 ,K— 1 ) ) 

UZET = U( I ,K) - U( I , K— 1 ) 

VZET = V( I , K ) - V( I , K-1 ) 

AZET = AA ( I , K ) - AA(I . K-1 ) 

XZET = X(I.K) - X(I .K-1) 

ZZET = Z( I ,K) - Z( I , K-1 ) 

YAC = XXI • ZZET - ZXI • XZET 

YAC = 1 . / YAC 

XIX = ZZET • YAC 

ZETAX- - ZXI • YAC 

XIZ = -XZET • YAC 

ZETAZ= XXI • YAC 

CNM = .5 • (CMU(I.K) + CMU( I ,K— 1 ) ) 

UX = UXI • XIX + UZET • ZETAX 

VX = VXI • XIX + VZET • ZETAX 

AX = AX I • XIX + AZET • ZETAX 

UZ = UXI • XIZ + UZET . ZETAZ 
VZ = VXI • XIZ + VZET • ZETAZ 

AZ = AX I • XIZ + AZET • ZETAZ 

TXX = -(-4. • UX + 2. • VZ) • CNM / 3. 

TXY = CNM • (UZ + VX) 

TYY = -CNM / 3. • (-4. • VZ + 2. • UX) 

R4 = ( (U( I , K )+U( I , K— 1 ) ) *TXX+(V( I , K )+V( I , K— 1 ) ) *TXY) *0 . 5 
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60 



1 + CNM / PR/ (GAL44A - 1 ) • AX 

S4 - ((U( I ,K)+U( I .K-1 )).TXY+(V( I .K)+V(I .K-1 ) ) .TYY) *0 5 
1 + CNM / PR / (GA>*4A - 1 . ) • AZ 

- 0 . 

- 0 . 

- 0. 

- (ZETAX • TXX + ZETAZ • TXY) 

= (ZETAX • TXY + ZETAZ • TYY) 

- (ZETAX • R4 + ZETAZ • S4) / 

2 . KMAXM1 

RH1 (K+1 ) - 
RH2(K+1 ) - 



70 



R4 
S4 

RH1 (K) 
RH2(K) 
RH3(K) 
RH4(K) 

DO 70 K 
DQ1 ( I .K) 
DQ2( I .K) 
DQ3( I .K) 
D04( I .K) 
CONTINUE 



/ YAC 
/ YAC 
YAC 



D0 1 ( I , K ) 
DQ2( I . K) 
D03( I , K ) 
DQ4( I , K ) 



RH3(K+1 ) - 

+ 1 ) * 



RH1 (K) 



RH4(K J 



RH2 (K J 
RH3(K) 
RH4(K) 



RETURN 

END 



SUBROUT I NE LOAD ( CPS . C L . CD . CM . A LFAS ) 

C0K440N/GR I D 1 /X ( 1 6 1 .41).Y(161.4l) 

COMy*ON/DGR I D/DT . I MAX , KMAX , I TEL. ITEU 
DIMENSION CPS(161) 

THIS SUBROUTINE COMPUTES THE INVISCID CONTRIBUTIONS 
TO LOADS ON THE AIRFOIL SURFACE 

IMAXM1 = IMAX - 1 
CL «= 0. 

CD « 0. 

CM = 0 . 

DO 400 I » ITEL . ITEU - 1 
XL = .5 • ( X ( I . 1 )+X( 1+1.1)) 

YL - .5 • (Y( I ,1 )+Y( 1+1,1)) 

DX «= X ( 1 + 1 . 1 ) - X ( I . 1 ) 

DY «= Y ( 1 + 1 • 1 ) ~ Y( I , 1 ) 

CPA - CPS(I+1) • .5 + CPS(l) • .5 
DCL = CPA • (-DX) 

DCD •= CPA • DY 
CL - CL + DCL 

CD «= CD + DCD 

400 CM = CM + DCD • YL - DCL • XL 

DCL - CL • COS(ALFAS) - CD • SIN(ALFAS) 

CD = CL • SIN(ALFAS) + CD • COS(ALFAS) 

CL = DCL 

RETURN 

END 



SUBROUT I N E WRAP ( 1 1 . J J . XS I NG . YS I NG . XP , YP . S0 . A0 . Y0 ) 

DIMENSION S0( 1 61 , 4) , Y0(41 , 4) , A0( 1 61 ,4),XP(1),YP(1) 

THIS SUBROUTINE UNWRAPS THE AIRFOIL AND STORES THE UNWRAPPED 
SURFACE IN ARRAYS A0 AND S0 . IT ALSO DETERMINES THE STRETCHING 
IN THE ETA DIRECTION. 

IMID - (II + 1) / 2 
DY - .8 / (JJ - 2) 

DO 1 J = 2 , JJ 
Y = FLOAT ( J-2) • DY 
1 Y0(J, 1) = 1 .25 • Y / (1 . - Y • Y) 

Y0( 1 , 1 ) = - Y0(3. 1) 

PI = 4. • ATAN ( 1 . ) 
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ANGL - PI + PI 
U = XP( 1 ) - XSING 

V « YP( 1 ) - YSING 
U - 1. 

V ® 0. 

I I M 1 - II - 1 
DO 2 I - 1 . II 
XII - XP( I) - XSING 
Y 1 1 *= YP ( I ) - YSING 

ANGL = ANGL + ATAN2( (U» Y1 1 -V.X 1 1 ) . (U*X1 1 +V.Y1 1 ) ) 

R - SORT (XI 1 **2 + Y1 1 »»2) 

U - XI 1 

V =Y1 1 

R — QAPT ( R ^ 

A0( I . 1 ) = R • COS( .5 • ANGL) 

2 S0( I . 1) - R • SIN( .5 • ANGL) 

C ! ! ! ! J I F OUTPUT OF UNWRAPPED COORDINATES IS DESIRED 
C WRITE (6. 1000) 

C WRITE (6.2000) ( I . A0( I , 1 ) . S0( I , 1 ) , I -1,11) 

RETURN 

1000 FORMAT( IX. 'UNWRAPPED COORDINATES IN THE TRANSFORMED PLANE') 
2000 FORMAT ( 1 5 , 2F20.8) 

END 

C 

c*.*.. 

c 

SUBROUT I NE TAB I NT ( XP , YP . XS ING . YS I NG . N) 

DIMENSION XP( 1 61 ) . YP(161 ) ,S0(161 ) ,A0( 161 ) 

C ! ! ! ! ! SMOOTH THE AIRFOIL SURFACE BY FINDING ADDITIONAL POINTS 
U = XP( 1 ) - XSING 

V = YP( 1 ) - YSING 
U = 1 . 

V * 0 . 

ANGL - 8. • AT AN( 1 . ) 

DO 1 I = 1 . N 

XII = XP( I ) - XSING 

Y 1 1 = YP(I) - YSING 

ANGL = ANGL + ATAN2( (U«Y1 1-V.X1 1 ) . (U«X1 1+V.Y1 1 ) ) 

R = SORT (XI 1 * *2 + Y1 1 •• 2) 

U = XII 

V = Y1 1 

R = SORT ( R ) 

A0( I ) = R . COS ( ANGL • .5) 

1 S0( I ) = R • SIN(ANGL • .5) 

DX =(A0(N)-A0( 1 ))/96. 

A0ST = A0( 1 ) 

DO 3 I = 1 .97 

XX = FLOAT ( I— 1 ) • DX + A0ST 

CALL TAINT(A0, S0.XX.YY.N.3.NER.MON) 

XP(I) = XX * XX - YY • YY + XSING 

3 YP( I ) - 2. . XX • YY + YSING 

RETURN 
END 
C 



c 

SUBROUT I NE T A I NT ( XTAB . FTAB . X . FX . N . K . NER . MON) 

DIMENSION XTAB( 1 ) , FTAB( 1 ) ,T(10) ,C( 10) 

C 

C NASA - AMES SUBROUTINE FOR POLYNOMIAL INTERPOLATION 
C OF A TABULATED FUNCTION 

C 

IF(N-K) 1,1.2 

1 NER = 2 
RETURN 

2 I F ( K— 9 ) 3,3,1 

3 I F (MON) 4.4.5 
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5 I F (MON-2) 6.7.4 
4 J - 0 

NM1 - N - 1 

DO 8 I * 1 . NM1 

IF(XTA8(I) - XTA8(I+1)) 9.11.10 

1 1 NER - 3 
RETURN 

9 J « J-1 
GO TO 8 
10 J = J+1 
8 CONTINUE 
MON = 1 

I F ( J ) 12 . 6 . 6 

12 MON « 2 

7 DO 13 I - 1 . N 

I F ( X - XTAB ( I ) ) 14.14.13 

14 J = I 

GO TO 18 

13 CONTINUE 
GO TO 15 

6 DO 16 I - 1 . N 

I F ( X-XTAB( I ) ) 16.17.17 

17 J ■= I 

GO TO 18 
16 CONTINUE 

15 J « N 

18 J * J — (K+1 ) / 2 
I F ( J ) 19,19,20 

19 J =■= 1 

20 M - J + K 

I F(M - N) 21 .21,22 

22 J = J - 1 
GO TO 20 

21 KP1 ■ K + 1 
JSAVE «= J 

26 DO 23 L - 1 . KP1 
C(L) = X - XTAB(J) 

T(L) - FTA8(J) 

23 J - J+1 

DO 24 J « 1 ,K 
I «= J+1 

25 T( I ) = (C(J).T(I)-C(I).T(J))/(C(J)-C(I)) 
I = 1 + 1 

I F ( I-KP1 ) 25,25.24 

24 CONTINUE 

FX - T(KP1 ) 

NER * 1 

RETURN 

END 



SUBROUTINE SING(N2,N.X.Z.XLE.YLE,TEA,TES.XSING. YSING.CHD) 



THIS SUBROUTINE COMPUTES SINGULAR POINT LOCATIONS. 

DIMENSION X(2) . Z ( 2) 

NU = N2 
N1 - N2 + 1 
N3 = N2 - 1 
XI = X(N1 ) 

Z1 = Z(N1) 

X2 = X(N2) 

Z2 = Z(N2) 

X3 = X(N3) 

Z3 - Z(N3) 

D1 - X2 •• 2 - XI •• 2 
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D2 = Z2 •• 2 - Z1 •• 2 

03 « X2 - XI 

04 « Z2 - Z1 

D5 - X3 •• 2 - XI • • 2 

06 - Z3 •• 2 - Z1 •• 2 

07 «= X3 - XI 

08 - Z3 - Z1 

B - (07 • ( D1 + 02) - 03* (05+06) )/(2. »(D7* 04-03*08)) 
I F( ABS(D3) . LT . ABS(D7) ) GO TO 10 



A “ 


(01 + 


02-2. 


• B 


• 04) 


/ 


(2. - 


► 03) 


GO 


TO 20 














10 A = 


(05 + 


06-2. 


• B 


• 08) 


/ 


( 2. 


• 07) 


20 CONTINUE 














R = 


SORT ( (X2— A) •• 


2 + 


(Z2-B) 


• 4 


■2) 





XLE = X ( NU) 

YLE = Z(NU) 

CHD - X(1) - X(NU) 

A2 = ( Z ( 2)-Z ( 1 ) ) /(X(2) - X(1)) 
A1 - ( Z ( N ) — Z ( N — 1 ) )/(X(N)— X(N— 1 )) 
TES «= .5 • ( A1 + A2 ) 

TEA = A2 - A1 

TEA = TEA • 57.29578 

XSING * (A+XLE) /2 . 

YSING = (B+YLE) / 2. 

RETURN 

END 



SUBROUTINE AIRF0L( II, JJ, IT, IE) 

COMwtON/GRID1/X( 161 . 41 ) . Z( 1 61 . 41 ) 

CCA440N/YSYM/ISYM 

DIMENSION S0(161 ,4) ,A0( 161 ,4) ,Y0(41 ,4) . XP( 1 61 ) . YP( 161 ) , 

1 E( 1 61 ) ,F(161).B0(49) 

DATA (B0( I) . 1-1 ,32)/1 . , 1 .0414, 1 .0836,1 .1270.1.1715.1.2175,1.2651. 
11.3145,1 .3659,1 .4199,1 .4755.1 .5349.1 .5973.1 .6636, 1 .7342.1 .8099, 
21.8914.1. 9799 . 2 . 0764 ,2.1829,2.3012,2.4341.2. 5653 . 2 . 7597 , 2 . 9646 , 
33.2106.3.5141,3.9019,4.4219,5.1687,6. 3632 . 8 . 6809/ 

! ! !!! COMPUTE THE COMPUTATIONAL GRID POINTS BASED ON INPUT AIRFOIL SHAPE 
DO 8 I - 1 .32 
8 A0( I , 1 ) = B0( I ) 

READ (5.1) 

1 FORMAT (IX) 

READ (5.2) FNU.FNL.ZSYM 

2 FORMAT (3F1 0 . 0) 

ISYM - 0 

I F(ZSYM.NE.0. ) ISYM - 1 



II > 


= 157 




JJ : 


= 41 




IT i 


= 31 




IE . 


= 127 




I IP1 


= II 


+ 1 


1 1 Ml 


= II 


- 1 


I IJJ 


- II 


• JJ 


IIJJ2 = II 


• ( JJ- 


ILE 


* (IT 


+ IE ) , 


I STP 


=* 0 




NN 


= 5 




NRF 


- 0 




NOTAPE - 1 




PI 


= 4. 


• ATAN( 1 


NU = 


FNU 




NL * 


FNL 




N = 


NU + 


NL - 1 


READ(5. 1) 
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READ (5.333) (XP( I ) . YP( I ) . I -= NL , N) 

333 FORMAT (2F 10.0) 

9994 FORMAT (F20. 8) 

L - N + 1 

I F ( 2SYM .GT. 0. ) GO TO 9995 
L = NL + 1 
READ(5 , 1 ) 

READ (5,333) (XP( L-I ) , YP ( L-I ) , 1=1 , NL) 

GO TO 9996 

9995 K1 - L 

DO 16 I - NL . N 
K = K 1 — I 
XP(K) = XP(I) 

YP(K) = - YP( I ) 

16 CONTINUE 

9996 SCALE = 1. /(XP( 1 )-XP(NL) ) 

XX = XP(NL) 

YY = YP(NL) 

DO 9997 I = 1 , N 

XP( I ) = XP(I) • SCALE 

9997 YP( I ) = YP ( I ) • SCALE 

CALL S I NG ( NU , N , XP . YP . XLE . Z LE . T EA . T ES . XS I NG . YS I NG . CHD ) 

CALL TABINT(XP.YP,XSING,YSING,N) 

NBODY = IE + 1 - IT 
DO 6791 1=1 . NBODY 

L = I - 1 
E(IT+L) = XP(I) 

6791 F(IT+L) = YP ( I ) 

IEP1 = IE + 1 
SLOPT = TES • .75 
DO 438 I = IEP1 . II 
II *= I +1 - IE 
E( I ) = A0(I1 . 1) 

DXI * 1 . / 48. 

D = 4. / 3. • (E( I ) - .25) 

F( I ) = F ( I E ) + SLOPT • ALOG(D) / D 
L = IIP1 - I 
E(L) = E( I ) 

438 F(L) = F ( I T ) + SLOPT . ALOG(D)/D 

C WRITE (6.439) 

439 FORMAT (2X.3H I . 1 9X . 1 HX , 1 9X , 1 HY) 

C WRITE (6.37) (I.E(I),F(I),I «= 1 . II) 

CALL WRAP ( I I . J J , XS ING , YSING . E . F , S0 , A0 , Y0) 

C ! ! ! ! ! MAP GRID BACK TO PHYSICAL PLANE AND SHIFT TO QUARTER CHORD 
DO 10 J = 2 , JJ 
DO 10 I * 1 . II 

X(l.J-l) = A0( I , 1 ) * • 2 - (S0(1 , 1 )+Y0( J , 1 ) ) *»2 
1 - 0.25 

10 Z(l.J-l) = 2. • A0( I , 1 ) . (S0(I,1)+Y0(J,1)) 

JJ = JJ - 1 
RETURN 

37 FORMAT ( I5.2F20.8) 

END 



SUBROUTINE CLUSTR(DS) 

COMMON /GR I D 1 /X ( 1 6 1 .41).Z(l61,41) 

COMMON/DGR I D/DT . IMAX.KMAX, I TEL. ITEU 
DIMENSION S ( 4 1 ) , XP (41), YP (41 ) , R ( 4 1 ) 

THIS SUBROUTINE CLUSTERS A GIVEN X.Z GRID SUCH THAT THE FIRST POINT IS AT 
THE USER-SPECIFIED DISTANCE DNM1N 
! ! ! ! .'COMPUTE THE OLD STRETCHING 
DO 100 I = 1 . I MAX 

S(1) = 0. 

XP( 1 ) = X ( I . 1 ) 
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YP(1) - Z(I .1) 

DO 10 K - 2 , KMAX 
XP(K) - X( I , K) 

YP(K) - Z( I , K ) 

10 S(K) - SORT ( ( XP(K )— XP(K— 1 ) ) • »2+( YP(K)-YP(K— 1 ) ) «*2) 
1+S(K-1 ) 

SUMDX - S (KMAX) 

CALL STRTCHf SUMDX , DS , FI , KMAX . FACTOR) 

C WRITE (6.200) I .FACTOR 

R(1) - 0. 

DR « DS 

DO 20 K « 2 . KMAX 
R(K) - R ( K— 1 ) + DR 
DR = DR • FACTOR 
20 CONTINUE 

RLAST - 1 . / R(KMAX) 

DO 30 K = 2 . KMAX 

R1 = R(K) • RLAST • S(KMAX) 

CM!! ! RED I STR I BUTE THE CONSTANT-ETA LINES 

CALL TA I NT ( S , XP . R 1 ,XP1 , KMAX . 3 , NER . MON) 

X( I .K) = XP1 

CALL TA I NT ( S . YP , R 1 ,YP1 , KMAX . 3 . NER . MON) 

Z( I , K) -= YP1 
30 CONTINUE 
100 CONTINUE 
C WRITE (6.115) 

DO 110 I = 1 . I MAX 
DX = X( I .2) - X( I . 1 ) 

DY = Z( I . 2) - Z(I.I) 

DN = SORT ( DX»DX+DY »DY ) 

C WR I TE( 6 , 120) I . DX . DY . DN 

110 CONTINUE 
RETURN 

115 FORMAT ( 5X.6HNORMAL, IX, 8HD I STANCE. 3H AT . 4H THE.5H WALL./ 
1 ,5H I . 8X . 2HDX , 8X , 2HDY , 8X , 2HDN,//) 

120 FORMAT ( I5.3F10.5) 

200 FORMAT ( 1 5 , FI 0 . 5) 

END 



SUBROUTINE STRTCH(SUMDX,DX1 ,F1 ,N1 ,R) 

THIS SUBROUTINE DETERMINES A GEOMETRIC 
PROGRESSION OF GRID SPACING BETWEEN 1 AND N1 SUCH THAT 
SUM8DX) EQUALS SUMDX. THE RATIO BETWEEN SUCCESSIVE 
SPACINGS IS R. 

N - N1 - 1 
R = 1 .5 
El = 1 . E— 04 
E2 = 1 . E— 04 
DO 10 L * 1, 50 

F= ( R-1 ) • SUMDX - DX1 • (R**N— 1 ) 

FP = SUMDX - DX 1 • FLOAT (N) • R •• (N-1) 

RITER = R - F/ FP 

I F ( 1 . E-02 . LT . R I TER . AND . RITER . LT . 1 . ) RITER « 1. 

I F ( 1 . .LT.RITER.AND.RITER.LT. 100.) RITER=.01 
I F ( ABS( R-R ITER) . LT . R*E1) GO TO 1 
R - RITER 
10 CONTINUE 
R = 1 .0001 

DX1 = DZTOT/FLOAT (N1— 1 ) 

RETURN 
1 R= RITER 
RETURN 
END 
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SUBROUTINE EDDY(DALFA) 

COMMON/ FLOW/Q1 ( 161 .41 ) ,02( 161 ,41 ) ,Q3( 161 .41 ) ,Q4( 161 .41 ) 

COMMON /MUTUR/CMU (161.41) 

C0M-40N/SK I NC F/C F ( 1 6 1 ) 

COMMON/DGRID/DT . I MAX . KMAX , I TEL , ITEU 

COMMON/PAR/GAMvlA . REYREF . ALFA . ALFA1 . REDFRE . AMINF . ALFA I 
COMAON/GR I D 1 /X (161 , 4 1 ) , Z ( 1 6 1 .41) 

DIMENSION T I N ( 4 1 ) , TOUT ( 4 1 ) , Y( 41 ) 

C 

C INITIALIZE viscosity everywhere 

FACT 1 = DT « AMINF / REYREF 
CMUMAX = 100. • FACT 1 / DT 

DO 1 K = 1 . KMAX 
DO 1 1 = 1 . I MAX 
1 CMU( I , K) = FACT 1 

C THIS SUBROUTINE COMPUTES THE EDDY VISCOSITY BASED ON THE 

C BALDWIN-LOMAX TWO LAYER MODEL 

C 

DO 100 I = 2 . IMAX - 1 
UDIF - 0. 

FMAX =0.1 E-06 
YMAX = . 1 E-06 
FYMAX = 0. 

Y(1) = 0. 

UWALL = 0. 

IF( I ,GT. ITEU. OR. I . LE. ITEL)UWALL = SORT (02(1,1) • *2+Q3( I , 1 ) »»2)/ 
101(1,1) 

C COMPUTE TAU AT THE WALL 

UET = 1 *(02( I ,2)/01 ( I .2) - 02(1.1 )/Q1 ( I . 1 ) ) 

VET = 1 . • ( 03 ( I , 2)/01 (1.2) - Q3( I . 1 )/Q1 (1,1)) 

XXI = X( 1 + 1 , 1 ) - X ( I — 1 . 1 ) 

ZXI = Z( 1+1 . 1 ) - Z( 1-1 . 1 ) 

XET = 4. • X ( I . 2 ) - 3. • X ( I , 1 ) - X( I , 3) 

ZET = 4. • Z( 1 , 2) - 3. • Z ( I , 1 ) - Z( I . 3) 

XXI = .5 • XXI 
ZXI = .5 • ZXI 
XET = .5 • XET 
ZET = .5 • ZET 

YAC = 1. / (XXI • ZET - ZXI • XET) 

OMEGA = (UET • XXI - VET • ZXI ) • YAC 

TWALL = AMINF • OMEGA / REYREF 
CF( I ) = 2. • TWALL / (AMINF. *2) 

FACT = SORT (01(1,1) • ABS( TWALL) ) .REYREF/ (26 . * AM INF) 

DO 10 K = 2 , KMAX— 1 

UXI = (Q2( 1+1 ,K)/Q1 ( 1+1 ,K) - 02(1-1 ,K)/Q1(I-1 ,K)) 

VXI = (Q3( 1+1 ,K)/Q1 ( 1+1 ,K)-03( 1-1 ,K)/Q1 ( 1-1 ,K)) 

UET = (02( I . K+1 )/01 ( I , K+1 )-02( I . K— 1 )/Ql ( I ,K-1 ) ) 

VET = (03(1 , K+1 )/Q1 ( I ,K+1 )-Q3( I . K— 1 )/Q1 ( I . K— 1 ) ) 

XXI = X( 1+1 , K) - X( 1-1 , K ) 

ZXI = Z( 1+1 .K) - Z( 1-1 , K ) 

XET = X( I , K+1 ) - X( I , K-1 ) 

ZET = Z( I , K+1 ) - Z ( I , K- 1 ) 

YAC = 1. / (XXI . ZET - ZXI • XET) 

OMEGA = A8S(UET.XXI+VET.ZXI-UXI*XET-VXI.ZET) • YAC 
UDIF = SORT (Q2( I , K) • .2+03 (1,K).»2)/Q1(I,K) - UWALL 
IF(ABS(UDIF) .GT.UDIFMAX) UD I FMAX = ABS(UDIF) 

Y(K) = SORT ((X(I,K)-X(I,K-1))..2+(Z(I,K)-Z(I,K-1))..2)+Y(K-1) 

F = Y(K) * OMEGA 
IF((Y(K)*FACT) .GT.20. ) GO TO 31 

IF( 1 .GT. ITEL.AND. I . LT. ITEU) F = F • (1. - EXP(-Y(K) .FACT) ) 

31 CONTINUE 
C 

C MODIFIED TURBULENCE MODEL APPLY FOR SPECIFIED RANGE OF ANGLES WHERE 
C FY IS USED T 0 FIND THE SECOND PEAK VALUE OF F FUNCTION 

C 
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IF(AIFA.LT ALFAI . AND , DALFA GE . 0 . ) THEN 
FY = F « Y(K) 

I F( FY . GT . FYMAX) THEN 
FYUAX = FY 
FMAX = F 
YMAX “ Y(K) 

END IF 
END IF 

I F ( ALFA . GE .ALFAI . OR . DALFA . LT . 0 . ) THEN 
IF(F.GT.FMAX) THEN 
FMAX - F 
YMAX = Y(K) 

END IF 
END IF 

FCT = Y(K) • FACT 
I F ( FCT . GT . 20 . ) FCT = 20. 

FCT - ABS(FCT) 

EL = .4 . Y(K) • (1 . - EXP(-FCT)) 

TIN(K) - 01 ( I ,K) • EL • EL • OMEGA 
TIN(K) - ABS( TIN(K) ) 

10 CONTINUE 
KSWTCH - 0 
FWAKE = YMAX • FMAX 
FI = 0.25 • YMAX • UDIF ..2 / FMAX 
IF(F1 . LT. FWAKE) FWAKE = FI 
DO 20 K - 2 . KMAX - 1 
FKLEB = 0. 

IF(ABS(Y(K)/YMAX) . LT. 1 . E+04) THEN 

FKLEB - 1. / (1. + 5.5 • (0.3 • Y (K)/YMAX) •• 6) 

END IF 

TOUT (K) - .0168 • 1.6 • 01 (I, K) • FWAKE • FKLEB 
TOUT (K) = ABS(TOUT (K) ) 

I F (KSWTCH . NE . 0) GO TO 20 
IF(TIN(K) .GT. TOUT (K) ) KSWTCH - K - 1 
20 CONTINUE 

C ! ! ! ! ! TOTAL VISCOSITY IS SUM OF LAMINAR AND INNER/OUTER LAYER AS APPROPRIATE 
DO 30 K = 2 , KMAX - 1 
I F(K . LE. KSWTCH) THEN 

CMU(I.K) - DT . (AMINF/REYREF + ABS(TIN(K))) 

ELSE 

CMU(I.K) - DT • (AMINF / REYREF + ABS(TOUT(K) ) ) 

END IF 
30 CONTINUE 
100 CONTINUE 
RETURN 
END 



SUBROUTINE RESI (OMEGA) 

C0M40N/PERTR/D01 (161 .41),DQ2(161 ,41),DQ3(161 .41),DQ4(161 ,41) 
C0MM0N/GRID1/X(161 ,41 ) ,Z(161 ,41 ) 

C0M40N/DGR I D/DT , IMAX.KMAX, I TEL, ITEU 

C0M40N/FL0W/Q1 (161 , 41 ) , 02 ( 1 61 , 4 1 ) ,03(161 .41 ) ,04( 1 61 . 41 ) 
C0M40N/PAR/GAM4A. REYREF, ALFA. ALFAI . REDFRE . AMI NF , ALFA I 
DIMENSION RHS( 161 ,4) 

XTAU( I , K ) - OMEGA . Z(I ,K) 

YTAU(I.K) - - OMEGA • X ( I , K) 

THIS SUBROUTINE COMPUTES THE RESIDUAL ON THE RIGHT HAND 
SIDE ARISING FROM THE EULER- PART OF THE GOVERNING EQUATIONS 

FLUX TERMS WITHIN THE XI- DERIVATIVE 
DO 100 K = 2 , KMAX - 1 
DO 10 I = 1 , IMAX 

UCON = (Q2( I ,K)/Q1 ( I ,K) ) . ( Z ( I . K+1 )-Z( I . K-1 ) ) 

1 - (03 ( I . K )/01 ( I , K) ) • ( X ( I ,K+1 )— X( I , K— 1 ) ) 

UCON =0.25 • DT • UCON 
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x i t - - xtau(i.k) • ( z( i , k+ i )- z ( i , k-i ) ) 

1 + YTAU(I.K) • (X(I,K+1) - X( I ,K-1 ) ) 

XIT - XI T . DT • 0.25 

UCON - UCON + XIT 

RHS( 1,1) - UCON • 01 ( I ,K) 

R - 1 . / 01 ( I . K ) 

P = (GAA44A— 1 . ) . (04(1. K) - .5 • R *(Q2(I ,K).»2+ 

1 03(1. K). *2)) 

RHS( 1.2) - 02(1, K) • UCON + P . DT • 0.25 • (Z(I.K+1) - Z(I.K-I)) 

RHS( 1,3) - 03(1. K) • UCON - P • DT . 0.25 • ( X( I , K+1 )-X( I . K-1 ) ) 

RHS( 1,4) = UCON • (04(1, K)+P) - XIT • P 

10 CONTINUE 

DO 11 1-2 . I MAX - 1 

DQI(I.K) = DOI(I.K) - RHS (1+1,1) + RHS( 1-1,1) 

DQ2 ( I , K ) = D02(I.K) - RHS(I+1,2) + RHS(I-1,2) 

DQ3 ( I , K ) = D03 ( I . K ) - RHS(I+1.3) + RHS(I-1.3) 

11 DQ4(I.K) - D04 ( I , K ) - RHS (1+1 .4) + RHS(I-1,4) 

100 CONTINUE 

FLUX TERMS WITHIN THE ETA- DERIVATIVE 

DO 200 1=2. IMAX - 1 
DO 20 K = 1 . KMAX 

VCON = (02( I ,K)/01 (I ,K) ) • (Z( 1-1 ,K)-Z( 1+1 ,K)) 

1 +(Q3(I ,K)/Q1 (I ,K) ) • (X( 1+1 ,K)-X( 1-1 ,K)) 

VCON = VCON • 0.25 • DT 

ETAT - -XTAU( I ,K) • (Z( 1-1 ,K)-Z( 1+1 ,K) ) -YTAU(I.K)* 

1 (X( 1+1 , K )— X (1—1 , K ) ) 

ETAT = ETAT • 0.25 • DT 
VCON *= VCON + ETAT 
RHS(K . 1 ) = VCON • 01 ( I ,K) 

P = (GAAKA-1.) . (04(1, K) - 0.5 • (02 ( I , K ) • .2+03 ( I . K ) • • 2 )/Q1 ( I , K ) ) 
RHS(K .2) = VCON . 02(1, K) +P • DT . .25 • (Z( 1-1 ,K)-Z( 1+1 . K) ) 
RHS(K , 3) = VCON . 03(1, K) + P • DT • .25 • (X( 1+1 . K) - X(I-I.K)) 
RHS(K , 4 ) = VCON • (04(1 , K )+P) - ETAT • P 

20 CONTINUE 

DO 21 K = 2 . KMAX - 1 

DQI(I.K) - DQI(I.K) - RHS(K+1 , 1 ) + RHS(K-I.I) 

DQ2 ( I . K ) = D02 ( I . K ) - RHS(K+1,2) + RHS(K-1,2) 

DQ3( I , K) = D03 ( I , K ) - RHS(K+1,3) + RHS(K-1,3) 

21 DQ4( I . K) = D04 ( I , K ) - RHS(K+1 , 4) + RHS(K-1 ,4) 

200 CONTINUE 

RETURN 

END 



SUBROUT I NE ROTGR I D ( X . Z , I MAX , KMAX . DALFA ) 

ROTATE GRID IN THE CLOCKWISE DIRECTION BY AN AMOUNT DALFA 
DIMENSION X( 1 61 ,41) , Z ( 1 61 ,41 ) 

CA = COS (DALFA) 

SA = - SIN(DALFA) 

DO 20 K - 1 . KMAX 

DO 20 I - 1 , IMAX 

XI = X( I ,K) 

Z1 = Z( I . K ) 

X( I , K ) = XI • CA - Z1 • SA 

20 Z( I , K ) - Z1 • CA + XI • SA 

RETURN 
END 



SUBROUTINE CPPL0T(I1 , I 2 . FMACH , X , Y . CP) 

THIS SUBROUTINE PLOTS CP AT EQUAL INTERVALS IN THE MAPPED PLANE 
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C0M40N/SK I NC F/C F ( 1 6 1 ) 

DIMENSION KODE(4) . L I NE( 90 ) , X ( 1 61 ) . Y( 1 61 ) ,CP( 1 6 1 ) 
DIMENSION CFX(3,49) ,CFY(3,49) ,CPX(3.49) .CPY(3.49) 

DATA KODE/1H , 1 H+ . 1 HI . 1H*/ 

. WRITE (6,2) 

. 2 FORMAT (50H0PLOT OF CP AT EQUAL INTERVALS IN THE MAPPED 

. 1 1 0H0 X/C , 1 0H CPU , 1 0H CFU , 1 0H 

CP0 - (1. + .2 • FMACH **2) •• 3.5- 1. 

CP0 * CP0 / ( .7 . FMACH »»2 ) 

K0 ■ 30. • CP0 +4.5 
IMIN - (12-11 )/2 + II 
I LOW = 2 • IMIN 
I COUNT = 0 
CHD=X ( 1 1 ) - X( IMIN) 

DO 12 I - 1 .90 
12 LINE(I) «= KODE( 1 ) 

DO 34 I - IMIN , 12 

K = 30. • (CP0 - CP(I)) + 4.5 

K 1 = 30. • (CP0 — CP( I LOW— I ) ) + 4.5 

IF(K.LT.I) K ■= 1 

I F (K 1 . LT. 1 ) K 1 - 1 

IF(K.GT.90) K «» 90 

I F ( K 1 . GT . 90) K1 - 90 

L I NE( K0) «= KODE(3) 

LINE(K) - KODE( 2 ) 

LINE(K1 ) = KODE(4) 

XOC *= (X(I) - X( IMIN) ) / CHD 
• WRITE (6.610) XOC. CP(I).CF(I),CF(ILOW-I). LINE 

L I NE( K 1 ) - KODE( 1 ) 

34 LINE(K) ■= KODE( 1 ) 

C* • • GENERATE PLOT3D CP AND CF PLOTTING FILES 
DO 500 I-IMIN, 12 

XOC = (X(I) - X( IMIN) )/CHD 
I COUNT « I COUNT + 1 
CPY( 1 . I COUNT) - 0.000000 
« 0.000000 

- CP( I LOW - I) 

- CF( I LOW— 1 ) 

- CP(I) 

- CF( I ) 

-= XOC 

- XOC 
= XOC 

- XOC 

- XOC 

- XOC 



PLANE/ 
CFL ) 



500 



610 



/EOF 



CFY( 1 . ICOUNT) 
CPY(2 . ICOUNT) - 
CFY(2. ICOUNT) - 
CPY(3. ICOUNT) - 
CFY( 3 . ICOUNT) =■ 
CPX(1 .ICOUNT) - 
CFX( 1 . ICOUNT) - 
CPX( 2 . ICOUNT) = 
CFX(2 , ICOUNT) - 
CPX(3. ICOUNT) - 
CFX(3. ICOUNT) - 
CONTINUE 
I DM « 3 

JDM = 12 - IMIN + 1 
WR I TE( 50) I DM , JDM 
WRITE(50)((CPX(I ,J). 

+ ((CPY(I.J). 

WR I TEf 60) I DM . JDM 
WR I TE( 60) ( (CFX( I ,J). 

+ ((CFY(I.J). 

RETURN 

FORMAT (4F1 0 . 4 , 90 A 1 ) 
END 



•1 . I DM) . J. 
= 1 . 1 DM) . J« 

= 1 . IDM) . J. 
«1 . IDM) . J- 



1 .JDM) . 
1 .JDM) 

1 .JDM) , 
1 .JDM) 
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TITLE. 

NACA 0012 AIRFOIL 



I MAX: 


KMAX : 


OT: 


WW. 


ALFA: 


ALFA1 : 


ALFA! 


161 


41 


.005 


5. 


15.00 


10.00 


19.00 


I TEL : 


ITEU : 


REYREF 


: DNMJN. 


TSTART : 


FORMAT : 


RSTRT 


31 


127 


3.45 


.00005 


-1 .0 


3.0 


TRUE 


CSTP : 


CPLT : 


NSTP • 


PSTP: 








0 


0 


1000 


50 









FNU : FNL : FSYM: 

33. 33. 1. 

X: Y: 

0 . 0 . 

0.0050 .01153 

.0125 .01894 

.0250 .02615 

.0500 .03555 

.0750 .04200 

.1000 .04683 

.1500 .05345 

.2000 .05737 

.2500 .05941 

.2600 .05962 

.2700 .05978 

.2800 .05992 

.2900 .05999 

.3000 .06000 

.3100 .05999 

.3200 .05992 

.3300 .05980 

.3400 .05965 

.3500 .05947 

.4000 .05803 

.4500 .05581 

.5000 .05294 

.5500 04952 

.6000 .04563 

.6500 .04137 

.7000 .03664 

.7500 .03161 

.8000 .02623 

.8500 .02055 

.9000 .01448 

.9500 .00807 

1.0000 .00126 



REDFRE. AM ] NF 
0.151 .283 

PITCH: PLUNGE: 
TRUE FALSE 
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